[settylab-dotto-bot]

GPSettings.param_scheme chooses which cells the length scale is estimated from. In da() it also covers d and mu. param_scheme changes no default behaviour. The default, None, resolves to each entry point's current behaviour: "condition1" in de() and "separate" in da(). This PR also removes sync_parameters, which param_scheme="pooled" replaces. It also fixes #32, which knowingly changes da() z-scores and PTPs when n_landmarks is at least the number of cells (see below).

Fixes #32

kompot.de(adata, "cond", "A", "B", gp=GPSettings(param_scheme="symmetric"))
kompot.da(adata, "cond", "A", "B", gp=GPSettings(param_scheme="pooled"))

value	meaning
None	the entry point's default: "condition1" for de(), "separate" for da()
"condition1" / "condition2"	estimate from one condition's cells and share the value; asymmetric by design
"symmetric"	estimate per condition and share a size-weighted combination; swap-invariant
"pooled"	estimate from both conditions stacked together and share
"separate"	each condition estimates its own

An explicit value (ls=, or d= / mu= passed to da()) pins that parameter under any scheme.

Default behaviour: bit-identical to main except in the #32 regime

Base 19ee1a1 (the merge base on main) and head 812c15b were both exported with git archive, selected by PYTHONPATH (asserted via kompot.__file__), and run in one Slurm job. Every entry point was called at its defaults: de(), da(), da(gp=GPSettings()), DifferentialExpression().fit/predict, DifferentialAbundance().fit/predict, the deprecated compute_differential_*() wrappers, and the CLI (kompot de and kompot da with flags only, and with each tree's shipped config template). The DE template fails identically on both trees, a pre-existing defect: the DE CLI does not route its landmarks: null into GPSettings.

• At 812c15b (before the DA: uncertainty is numerically unstable when n_landmarks >= n_combined (the GPSettings default on <5000 cells) #32 fix): 75 of 75 arrays bit-identical. Determinism control, base against a second base run: 75 of 75.
• At 6606ee2: 70 of 75. The 5 that differ are all da(gp=GPSettings()) on 450 cells, which is in the DA: uncertainty is numerically unstable when n_landmarks >= n_combined (the GPSettings default on <5000 cells) #32 regime: z up to 6.3e-2, PTP up to 3.7e-2, fold changes and densities up to 4.7e-6. Direction labels are unchanged.
• Positive controls differ from the base default: de() with "symmetric" on 8 of 9 arrays, da() with "pooled" on 5 of 6, and a relative 1e-7 change to ls_factor on 8 of 9 (max 2.1e-07).
• The only difference at the defaults is metadata: the run record in uns gains gp.param_scheme.

Why de() needs a scheme

Condition 2 inherits condition 1's fitted length scale, so the default result depends on argument order. On an exchangeable null (2,917 cells, 236 runs), exchanging the conditions moves the fraction of p < 0.05 by 0.06–0.08 at 2:1 and beyond under the default. Under "symmetric" or "pooled" it moves by at most 0.0061. "separate" roughly doubles the false-positive tail (0.322 vs 0.172 at 2:1) and is documented as diagnostic only. "pooled" shrinks the length scale purely with cell count. "symmetric" does not, because it uses the size-weighted geometric mean of the per-condition estimates. This PR does not change which scheme is the default.

de(X, Y, "condition1") and de(Y, X, "condition2") are the same computation with the labels exchanged, provided both runs use the same landmarks. Automatic landmarks are row-order dependent, so pass one landmarks array or n_landmarks=0. fit() warns when "condition2" runs with automatic landmarks.

Differential abundance

• At "separate" without landmarks, da() is swap-equivalent, and a test covers this.
• "condition1" / "condition2" reproduce that condition's own estimator bit for bit and share the result.
• "symmetric" pools within-condition statistics: a geometric-mean ls, a d averaged over each condition's query cells, and mu over both conditions' nearest-neighbour distances at that d. It is bit-identical under a swap by construction. Its calibration has not been benchmarked.
• "pooled" is order-dependent, as sync_parameters=True always was. Two remedies are needed, and neither substitutes for the other: pin d, mu and ls, and share one landmarks array. See DifferentialAbundance.fit.
• Fixed: da() silently ignored GPSettings.ls; it now honours it. It also warns when the expression-only sigma, eps or use_empirical_variance is set. Its regression test fails before this change and passes after it.

Fixed: #32, abundance uncertainty depended on argument order

At n_landmarks >= n_combined, DifferentialAbundance.fit used the stacked union [X_condition1; X_condition2] as the landmark array. mellon's uncertainty, a diagonal Laplace approximation in Cholesky-whitened coordinates, depends on the landmarks' row order, so da(X, Y) and da(Y, X) disagreed on z by a median of 82% (PTP 96%), while fold changes agreed to 1e-10. Each condition now gets its full GP on its own cells, the same as n_landmarks=None.

entry point, at its defaults	in the regime?	result
da() without gp; DifferentialAbundance(); compute_differential_abundance(); CLI flags-only or shipped template	no	bit-identical
da(gp=GPSettings(...)) on ≤ 5000 cells; CLI with any GP key but no n_landmarks	yes	z and PTP move to the full GP; orientations agree exactly
anything with fewer landmarks than cells; all of de()	no	unchanged

On 400 + 400 cells, z moves by a median of 34%, and fold changes and densities by at most 6.0e-4. The result matches gp=None to 3e-14. The regression test fails at 812c15b and passes here, and four mutants flip exactly as predicted.

The upstream part belongs to mellon. With fewer landmarks than cells, k-means centroid order still reaches the uncertainty: on the same fixture at 100 landmarks it differs by a median of 87%. Sharing one landmarks array removes the dependence, which is already documented and tested. An order-invariant uncertainty would have to come from mellon, and is proposed for consideration in settylab/Mellon#22.

Fixed: kompot dm reported a missing PCA only when Palantir was installed

run_dm imported Palantir before validating its input, so test_dm_missing_pca failed in any environment without that optional dependency. CI installs .[all] and passed. The import now comes after the input checks.

Removed

sync_parameters on DifferentialAbundance.fit and da(), with no alias. git grep -n sync_parameters finds only the CHANGELOG line recording the removal. A caller who still passes it gets TypeError: DensityEstimator.__init__() got an unexpected keyword argument 'sync_parameters'. param_scheme takes its argument position, so allow_single_condition_variance stays where it is on main.

Verification

• Full suite at da733b4: 1 failed, 2112 passed, 56 skipped. The only failure is the worktree-only provenance test (_provenance: sha does not resolve in a git worktree with a branch checked out #33), which does not fail in CI's checkout. test_dm_missing_pca has left the failing set.
• Mutants, with predictions written first: condition2 reading condition 1's cells, a non-commutative symmetric weight, and dropping gp.ls in da(). Each flips exactly the tests predicted.
• Environment: python 3.12.12, numpy 2.4.4, jax 0.10.0 (CPU), mellon 1.7.1, anndata 0.12.10.

CI

794bf7f bounds and diagnoses the intermittent 3.11 hang. It sets a 60-min job timeout and fail-fast: false, and adds pytest-timeout at 300 s with thread stack dumps, plus --durations=25. It does not cure the hang.

Not decided here

• Which scheme should be the default for each entry point.
• The n_landmarks default split in da(): None without gp, 5000 with it.

[settylab-dotto-bot]

Regression fixed — 487a02e

An independent skeptic review found a real defect in the first commit, and it would have shipped.

The ls_scheme gate excluded cov_func / cov_func_curry from the automatic path. Under the
default scheme that dropped the historical ls_for_model2 = self.model1.ls inheritance, so a
caller who supplied a custom kernel and never touched ls_scheme silently got "separate"
behaviour — the worst-measured scheme in the comparison this branch is based on. function_kwargs
is public and kompot/cli/de.py forwards into it, so the path is reachable from the shipped API.
It falsified this PR's claim that nothing changes unless ls_scheme is set.

The gate is now exactly the historical one, ls is None and "ls" not in function_kwargs. A
supplied kernel no longer disables scheme resolution — the resolved value reaches the custom kernel
through estimator_defaults["ls"], which is what the old inheritance did as well, so both
conditions keep a shared scale.

Verified against 19ee1a1 in a clean worktree, default scheme, ls=None:

function_kwargs	old model1.ls / model2.ls	new
none	92.589588 / 92.589588	92.589588 / 92.589588
cov_func_curry=Matern52	92.589588 / 92.589588	92.589588 / 92.589588

Three tests added for the curry path (default, shared schemes, separate); 20/20 green. The
reason 17/17 missed it is that no test exercised function_kwargs at all — only explicit ls.

The docstring's "ignored when ... cov_func is supplied" was true of the scheme and false of the
behaviour; it now says the resolved value is handed to the custom kernel.

Full skeptic findings are on settylab/dotto-nexus#310. Everything else in the review reproduced or
strengthened the measurements; this was the one merge-blocking item.

[katosh]

@settylab-dotto-bot also implement the "condition2" ls_scheme.

[settylab-dotto-bot]

Picked up — dispatched to a worker on this branch.

The shape it is building to, so the intent is on the record before the code lands: "condition2" is the mirror of "condition1", fitting the shared length scale from condition 2's cells and giving it to model 1. The reason it earns a place rather than being a fifth knob is that it makes the default's swap-dependence demonstrable in a single run — de(X, Y, ls_scheme="condition1") and de(Y, X, ls_scheme="condition2") should be the same computation with the labels exchanged, which is an exact equivalence a test can assert rather than a property a benchmark has to estimate.

That equivalence is also the reason the 236-run null benchmark is not being re-run. If the mirror property holds exactly, the published condition1 row of the swap-invariance table already describes condition2 under label exchange. The worker has been told to establish that by test first and to distrust it until then — it is precisely the convenient claim that rescues it from an expensive measurement, and the pipeline may not be symmetric under label exchange everywhere (fold-change sign, anything indexed by position rather than by label). Whatever does not mirror gets reported at the same prominence as what does, and if the property cannot be established the row will be marked unmeasured rather than derived.

It has also been pointed at the gate that 487a02e fixed, since a new scheme value touches it and that defect — a custom cov_func silently turning the default into separate — is exactly the kind that re-instantiates when the enumeration grows.

[settylab-dotto-bot]

"condition2" implemented — 3cfc904

The mirror is real and it is bit-exact, so the derivation it licenses is on the record with a
test as its warrant rather than an argument.

The equivalence. With identical landmarks, de(X, Y, ls_scheme="condition1") and
de(Y, X, ls_scheme="condition2") agree at max_rel = 0.000e+00 on condition1_smoothed /
condition2_smoothed (roles swapped), both posterior standard deviations,
mahalanobis_distances and neg_log10_ptp. fold_change, mean_log_fold_change and
fold_change_zscores are an exact sign flip — that is what fold_change = condition2 - condition1 means, and it is reported as a finding rather than tolerated by weakening the test.
Control numbers, so the sign-flip rows are a statement and not two arrays of zeros agreeing: the
unflipped comparison gives max_rel = 2.000e+00 (the signature of u = -v) at
max|fold_change| = 0.5054.

What does not mirror, and it is not about ls_scheme. fit() derives automatic landmarks
from np.vstack([X_condition1, X_condition2]), whose row order differs between the two
orientations. compute_landmarks is order-sensitive, and at random_state=None — the
DifferentialExpression default — it is not even run-to-run stable at a fixed order. With
automatic landmarks that is the only residual disagreement between the orientations (max_rel
1.5e-05 on the smoothed surfaces, 1.8e-06 on the Mahalanobis distances, against 0.000e+00 when
one landmark array is shared). A pre-existing, scheme-independent swap asymmetry that no
ls_scheme removes, and one this PR had not named. It is why the equivalence is stated as given
identical landmarks, and it deserves its own issue.

The swap-invariance row is derived, not re-measured. Tail mass is sign-blind, so
T(de(A,B,"condition2")) = T(de(B,A,"condition1")) and vice versa; swapping the two terms of an
absolute difference changes nothing. condition2's row therefore is condition1's row,
identically. The 236-run benchmark was not re-run. What is not derivable and is not claimed:
condition2's absolute tail mass per imbalance ratio — those are the default's numbers on the
swapped split, published only in aggregate.

The gate from 487a02e was checked, not assumed. It is ls is None and "ls" not in function_kwargs — value-blind, so condition2 inherits the fix rather than needing its own. That
is now a measurement (test_condition2_still_shares_when_a_cov_func_curry_is_supplied), and the
reason it is safe is structural: the scheme value is consumed one level down in arms that are
literal equality over disjoint strings, where adding a member cannot shadow an existing arm.

Two readers of ls_scheme that could not see it, both live on this branch before this commit:

1. params_get could not read it back — build_params_dict stores it nested at
   params["gp"]["ls_scheme"], and _LEGACY_MAP did not list it, so it resolved to None and
   params_match was False for every rerun, including a byte-identical one. The
   "participates in run-parameter matching" claim was half-true: it joined the list, and the list
   could not read it.
2. The DE CLI silently misrouted a flat ls_scheme config key. An unrecognised key is not
   rejected — it falls through to de(**function_kwargs) and is handed to mellon, so the
   documented knob was unreachable from a config file. gp_keys is now module-level
   GP_CONFIG_KEYS, asserted against the dataclass's own fields.

Tests and totals. tests/test_ls_scheme.py 20 → 33. Negative control: planting
_auto_ls(X_condition2, ...) → _auto_ls(X_condition1, ...) flips exactly the 7 mirror tests,
and the two condition2 tests that stay green are the two that should. Full suite 1 failed,
2061 passed, 56 skipped against 1 failed, 2048 passed, 56 skipped at 487a02e in a clean
worktree — same failing-ID set, not just the same count, and the single failure is an absent
palantir in this environment.

One implementation note worth a reviewer's eye: "condition1" resolves lazily (model 1 estimates,
model 2 inherits) and "condition2" cannot, since model 1 is fitted first. It resolves eagerly
through _auto_ls, so the mirror holds only while the eager and lazy estimates agree exactly.
They do, to the last bit — pinned by a test rather than assumed.

Tracking: settylab/dotto-nexus#310.

[settylab-dotto-bot]

kompot PR #22 — the "condition2" ls_scheme

Implemented the "condition2" ls_scheme on kompot PR #22 and established the mirror property by measurement: given identical landmarks, de(X, Y, "condition1") and de(Y, X, "condition2")...

Full report: https://github.com/settylab/dotto-nexus/blob/2da4755a4837d325e05498519ebcb9dae5a4b4eb/assets/22/kompot-310-ls_2026-09-21_165401_kompot-pr22-ls-scheme-condition2.md

[settylab-dotto-bot]

Corrections — 082a702, 55d4ef2, 2bff8e3

An independent skeptic pass found three defects in what I claimed, none in what was built.
Recording them plainly, because the earlier comment and the PR body both overstated the result.

1. The equivalence was stated unconditionally in three places, and the clause that makes it true
lived in a test name. settings.py, the fit() docstring and the CHANGELOG all asserted
de(X,Y,"condition1") == de(Y,X,"condition2") with no precondition; the CHANGELOG said "measured
bit-identical". Given identical landmarks appeared only in the PR body, a commit message, and
the name of a test. A boundary that lives in a test name is invisible to everyone it protects.

And the configuration that violates it is the library default: n_landmarks=5000, landmarks=None. Automatic landmarks come from the two conditions' cells stacked in the order
given, so the two orientations select different sets — silently, with no error and no log line.

Fixed in all three places, each now saying what happens when the precondition does not hold.
Plus enforcement, because documentation is not a mechanism: fit() logs a warning when
ls_scheme="condition2" runs with automatic landmarks, scoped to that scheme because the
cross-orientation comparison is its only purpose.

2. "It is small" was one configuration generalised. Landmark-fraction sweep, 800 cells,
RMS-relative gap ‖u−v‖₂ / ‖u‖₂ — ~45× on the smoothed surfaces and ~25× on the Mahalanobis
distances over a 10× fall in the fraction, with a shared-landmark control at 0.000e+00 at both
ends, which pins the growth on landmark selection rather than on the fraction. A default
5,000-landmark run on a large dataset sits at the low-fraction end — the worse end. The table is
in the PR body; the sentence is now scoped to what was measured.

3. A miss count is not a number, and my "7, 8, 7" was three samples dressed as a range.
Twelve independent draws at the default random_state=None:

comparison	draws missed
condition1_smoothed ↔ condition2_smoothed	12/12
condition2_smoothed ↔ condition1_smoothed	12/12
condition1_std ↔ condition2_std	8/12
condition2_std ↔ condition1_std	2/12
fold_change vs negated	12/12
mean_log_fold_change vs negated	1/12
fold_change_zscores vs negated	12/12
mahalanobis_distances	12/12
neg_log10_ptp	12/12

Per-draw totals 6 (×4), 7 (×6), 8 (×1), 9 (×1) — so the observed range is 6–9, and my
published 7–8 was wrong at both ends. Six occurs at the default, and it is the same six
comparisons random_state=0 gives, so a reader who runs once at the default and sees 6 learns
nothing about the seed either way. Six is a floor, not a typical value.

The docs now publish the localisation rather than a count: six always-miss, the two posterior
standard-deviation comparisons carrying nearly all the variance, and mean_log_fold_change
nearly exempt because averaging over cells cancels most of the landmark noise. Nearly — it
missed once in twelve, so "never" would have been the same mistake one level down.

And the presentation itself was part of the defect. The two comparisons carrying all the
instability are exactly the pair the PR table folded into one "condition1_std /
condition2_std, role-swapped" row. That fold made a reviewer enumerate eight comparisons where
the population is nine, dropping a passing one and producing a short count. The table is now
unfolded and numbered. A collapsed table row is not a population, and here the rendering hid the
member and the variance together.

A note on the new tests, because the failure mode generalises. The positive "it warns" test
failed on its first run: kompot/__init__.py sets propagate=False on the logger, so pytest's
caplog sees an empty record list. That does not fail loudly — it makes every "did not warn"
assertion pass for a reason unrelated to the code, and my three negative tests were vacuous
until the positive one flushed it out. Replaced with a handler attached directly to the kompot
logger, plus a positive control for the helper. Forcing the warning to fire unconditionally flips
exactly the three negatives.

That is worth a look beyond this PR: tests/test_groupwise_fdr_integration.py::test_no_length_mismatch_warnings
asserts an absence through caplog.records on kompot's logger and therefore cannot fail —
measured, caplog.records = 0 with a canary WARNING emitted immediately before. Not touched here,
because once that assertion can see records it may start failing, and that belongs in its own
change rather than inside a PR about length scales. tests/test_volcano_de_rendering.py already
carries a correct helper for this with a docstring explaining the trap.

Numbers throughout are my own runs with a stated metric, not transcriptions: the skeptic's counts
came from its clone and its own metric, and where we differ — its ~100× growth factor, its stable
7, its 0/10 for mean_log_fold_change — the figures published here are the ones measured with
the fraction varied alone and the draw repeated twelve times.

Unchanged: the implementation survived the attack. The mirror reproduces bit-identically in an
independent clone, the bit-identity is structural rather than platform luck, the cov_func gate
is not re-opened, and the derived swap-invariance row remains correctly labelled derived with its
warrant named. tests/test_ls_scheme.py 33 → 38.

Tracking: settylab/dotto-nexus#310.

[katosh]

@settylab-dotto-bot make sure differential abundance works equivalently. Also keep the changelog concise, similar to how we wrote it for prior releases.

[settylab-dotto-bot]

Both taken, and the abundance one is a real gap rather than a restatement — flagging that rather than just acknowledging.

Differential abundance. The condition2 mirror work and tests/test_ls_scheme.py establish equivalence for differential expression: de(X, Y, "condition1") == de(Y, X, "condition2"), signed quantities flip, unsigned statistics hold. Abundance is a separate code path and nothing so far shows it mirrors. It is being established to the same standard rather than assumed to inherit — result-level assertions, and the positive controls before the flip check (|fold_change|.max() > 1e-2 and not allclose), which is what keeps the test from passing vacuously on an all-zero result. If abundance turns out not to mirror, that is the finding and it will be reported as one rather than patched toward the expected answer.

Changelog. Trimming to match the register of the prior entries rather than to a guess at "concise". The landmark precondition currently in the entry is a real caveat — automatic landmarks are computed from the two conditions' cells stacked in the given order, so the two orientations select different landmark sets and the equivalence is only approximate under the defaults — so it moves to the docstring, where a user meets it at the point of use, rather than being dropped. The changelog entry will point at it.

One note on scope: the exact equivalence measured so far is conditional on identical landmarks in both runs. That condition is not satisfied by the defaults, and the same question applies to abundance.

[settylab-dotto-bot]

Differential abundance: already equivalent — c907da9

make sure differential abundance works equivalently. Also keep the changelog concise

Abundance already works equivalently, and it needs no ls_scheme analogue. That is the
answer, and it is a better one than adding a knob would have been.

DifferentialAbundance.fit builds both density estimators from the same estimator_defaults
dict, and there is no ls_for_model2 = self.model1.ls analogue anywhere in the module — no
condition-1 inheritance to be asymmetric about. The defect ls_scheme exists to expose in
differential expression has no counterpart in abundance. (A claim about the code as read at
082a702.)

The whole argument rests on that one line, so it was measured, not asserted: **d unpacks into a
fresh dict, but the values are shared references, and an in-place mutation of the landmarks array
by estimator 1's fit would make the result luck rather than structure. estimator_defaults is
byte-identical after constructing estimator 1 and again after fitting it; landmark object
identity stable; both estimators' landmarks equal the dict's.

Measured at the abundance defaults (n_landmarks=None, sync_parameters=False, 400 + 400 cells,
2-D, seed 0) — the configuration named because it is the reason, since with no landmarks there is
no np.vstack row-order dependence to be sensitive to:

quantity	comparison	max_rel
log_density_condition1 / condition2	role-exchanged	0.000e+00
log_fold_change_uncertainty	vs swapped	0.000e+00
neg_log10_fold_change_ptp	vs swapped	0.000e+00
log_fold_change, log_fold_change_zscore	vs negated swapped	0.000e+00
log_fold_change_direction	vs label-exchanged swapped	800/800

One dataset, one configuration. It generalises through the mechanism, not through the zero.

Two configurations are only approximate, and both are relatively worse than expression's
1.5e-05 landmark gap — sync_parameters=True at 1.22e-02 on the log fold change, and
automatic n_landmarks=40 at 5.6e-02 on the uncertainty. Same root cause with a wider blast
radius: np.vstack([X_condition1, X_condition2]) is row-order dependent, and under
sync_parameters it feeds compute_d_factal / compute_mu / compute_ls as well as
compute_landmarks. So the landmark precondition documented for expression is in scope for
abundance and has a second instance there. Sharing one landmarks array restores exactness,
and that is tested.

The control is why these tests are worth anything

My first two direction checks read 160/160 and 600/600 agree. Both were vacuous — every
cell was 'neutral' (max|lfc| 0.07, then 0.42, against the package's own default threshold of
1.0), and two all-'neutral' arrays agree perfectly whether or not the labels are swapped. A
single arm reading 600/600 looks like a triumph; only the unswapped arm exposed it.

Why a ported test could not have caught it, which is the part that generalises: differential
expression has no categorical output, so a control built for the expression family structurally
under-covers the member being added. A ported shape inherits the green without earning it. The
direction test now asserts both arms and requires 'up' and 'down' to occur; separation was
raised until |log_fold_change| clears the default threshold, and no threshold was tuned.

Changelog

Cut from 71 lines to 42, with the register derived from the 0.8.0 entries in the file rather than
guessed — bolded lead, one to three sentences. The landmark precondition is not deleted: it lives
in DifferentialExpression.fit (a .. warning:: block) and GPSettings.ls_scheme, where a user
meets it at the point of use, and the changelog points there.

Also, from the wider sweep

tests/test_groupwise_fdr_integration.py::test_no_length_mismatch_warnings asserts
len(length_warnings) == 0 over caplog.records, and cannot fail: kompot/__init__.py sets
propagate=False, so caplog sees an empty list. Measured — caplog.records = 0 with a canary
WARNING emitted immediately before. An AST sweep of tests/ finds 12 functions taking caplog, 3
with absence-shaped assertions, 2 of which are unrelated false positives; that one is the only real
instance. Note its assertion text contains no caplog — the binding is two statements upstream,
so only data flow finds it.

Not fixed here: once it can see records it may start failing, which would be a live disclosure in
the groupwise-FDR path landing inside a PR about length scales. It wants its own change.
tests/test_volcano_de_rendering.py already carries a correct helper for this with a docstring
explaining the trap.

One measurement outstanding, named rather than glossed: the last full-suite run predates
tests/test_da_swap_equivalence.py, so its 8 tests pass in isolation but not yet inside a
whole-suite run. No library source changed in the commits that added them, but a new test file can
affect others through import or fixture side effects, so the re-run is pending rather than assumed
unnecessary.

Tracking: settylab/dotto-nexus#310.

[settylab-dotto-bot]

Skeptic pass — differential-abundance swap equivalence, kompot PR #22 @ c907da9

Verdict: check — the differential-abundance swap equivalence holds and is wider than measured, but two published sentences over-reach and a third, pre-existing numerical issue surfaced.

Full report: https://github.com/settylab/dotto-nexus/blob/187bfb922a74eb73d86fcccdadca0eb044de1d35/assets/22/kompot-310-dask_2026-09-22_101630_kompot-pr22-da-swap-skeptic.md

[settylab-dotto-bot]

Retraction — two over-reaching claims corrected in f72573e

An independent review pass found two claims on this PR that were wrong, both mine, and both the
same move: a correct local measurement generalised past the condition that made it true.
Retracting them here because one of them is in a comment above, which cannot be edited without
the edit being invisible.

---

1. "Passing one landmarks array to both runs restores exactness" is FALSE for the
sync_parameters=True half.

I wrote that in the PR body, in my 16:41:56Z comment, in CHANGELOG.md, and in the test module
docstring. It is measured for the landmark half only. It is now measured false for the
other: sync_parameters=True with 40 shared landmarks reads 4.784e-03, against a 4.720e-03
reference with no landmarks at all. Indistinguishable — and the reference arm is what makes that a
number rather than an impression. mu and ls come from the union's nearest-neighbour distances,
which landmarks are not an input to.

The correct remedy, with its size boundary in the same sentence: pass mu and ls
explicitly — and d as well when the combined data exceeds 500 cells, because
compute_d_factal subsamples 500 indices above that and becomes row-order dependent itself.
Below 500, mu and ls suffice. Passing all three is always correct.

configuration	combined cells	log_fold_change
d + mu + ls explicit	800	0.000e+00
mu + ls only	800	6.085e-04 (d reads 1.4972 vs 1.4983)
mu + ls only	400	0.000e+00

That boundary was one arm from being missed, and the near-miss is the same defect again. My
own earlier note that "d is identical across orientations" was measured at 160 combined cells.
True there; false above 500. Had the mu+ls arm been run at a single size, "pass mu and ls"
would have shipped as the remedy and been wrong for nearly every real dataset — a corrected
prescription that was itself unscoped.

---

2. "Differential abundance is swap-equivalent at its defaults" did not say WHICH defaults, and
the entry point users configure through does not have them.

da() takes n_landmarks=None only when gp is None
(anndata/differential_abundance.py:104), and any GPSettings carries n_landmarks=5000
(settings.py:82). So da(adata, gp=GPSettings(...)) lands in the automatic-landmark branch that
the same paragraph calls approximate. GPSettings being the documented way to configure da(),
the configuration that is exact is arguably the one a user is least likely to be in. The
measurement was never wrong; the sentence around it claimed a reach it did not have.

And at that configuration it is not merely "approximate": when n_landmarks >= n_combined the
landmark set becomes the whole union and the uncertainty-derived quantities are numerically
unstable between orientations. That is pre-existing, traces to the 0.8.0 gp_type="fixed"
change, and is not caused by this PR — tracked separately as #32.

---

What is unchanged

The measurements themselves. The abundance mirror is still bit-exact at
n_landmarks=None, sync_parameters=False; condition2 still mirrors condition1 bit-exactly
under shared landmarks; the swap-invariance row is still derived with its warrant named. Nothing
retracted here is a number — all three are sentences that carried a number past its scope.

Three instances of one habit on one PR — "at its defaults", "shared landmarks restore
exactness", and "d is identical" — which is why they are being reported together rather than
patched quietly. The pattern is worth more than any of them individually: a measurement is true
of the configuration it was taken at, and the sentence around it is where that gets lost.

Corrected in f72573e: CHANGELOG.md, the PR body above, tests/test_da_swap_equivalence.py
(module docstring and test_shared_landmarks_restore_the_exact_swap, now scoped to the landmark
half). The test file is otherwise unchanged — verified AST-identical with docstrings stripped.

Tracking: settylab/dotto-nexus#310.

[katosh]

@settylab-dotto-bot this whole setup makes the API of da and de somewhat inconsistent, both allow parameter synchronization but under different mechanisms. What would be the cleanest way to make this consistent, easy to understand and well integrated into our API?

[settylab-dotto-bot]

Making de and da parameter sharing consistent: one field in GPSettings

this whole setup makes the API of da and de somewhat inconsistent … What would be the cleanest way to make this consistent, easy to understand and well integrated into our API?

Recommendation: a single GPSettings field, read by both de() and da(), that says where every hyperparameter the model estimates from cells comes from. It replaces both ls_scheme and sync_parameters. No new concept, no second axis:

@dataclass
class GPSettings:
    ls: Optional[float] = None
    ls_factor: float = 10.0
    param_scheme: Optional[str] = None   # "separate" | "condition1" | "condition2" | "symmetric" | "pooled"
    ...                                  # None = this entry point's default (the module's existing convention)

kompot.de(adata, "cond", "A", "B", gp=GPSettings(param_scheme="symmetric"))   # shares ls
kompot.da(adata, "cond", "A", "B", gp=GPSettings(param_scheme="pooled"))      # shares d, mu, ls

What the scheme covers is decided by the model rather than by the user. In de it covers ls, the only hyperparameter FunctionEstimator estimates from cells: it has no d, and its mu is a constant prior mean of 0. In da it covers d, mu and ls. An explicit value (ls=, or d=/mu= for da) pins that parameter and takes precedence, as it already does in both code paths.

Every existing value maps exactly, so no numbers move:

today	becomes
de: ls_scheme=X, any of the five	param_scheme=X (the same code, renamed)
de: default ("condition1")	param_scheme=None → "condition1"
da: sync_parameters=False (default)	param_scheme=None → "separate"
da: sync_parameters=True	param_scheme="pooled"

The two trade-offs that decide it:

1. One concept instead of two, and what you give up is "share ls but not d/mu" in da. I weighed the alternative of two fields, a source plus a scope (none|ls|all). In de the scope field would have only one effective value, because ls is all there is. In da nobody has measured whether sharing ls alone is ever what you want. A scope field can be added later without breaking anything. Removing one after release cannot be done so cleanly.
2. Defaults stay different, but now in plain sight. Under one vocabulary the real inconsistency is one line of documentation: de defaults to "condition1" and da defaults to "separate". Whether da should share by default is a scientific question. The harm from unshared length scales is measured for DE only, so I would not settle it as a side effect of tidying the API. A None default is what lets the field keep each entry point's current behaviour. It is also the convention the settings.py module docstring already states: "None means 'use kompot's default'" (kompot/settings.py:10).

One timing point: ls_scheme is still under [Unreleased], so renaming it now costs nothing. After the next release it needs a deprecation cycle. sync_parameters has shipped, so it becomes a deprecated alias on DifferentialAbundance.fit that maps True/False onto the table above.

Evidence (code at f72573e; probe on synthetic data: 700 + 300 cells, 3-D, seed 0)

• sync_parameters=True computes the same length scale as pooled, to the bit. da uses compute_ls(compute_nn_distances(vstack(...))) * ls_factor (kompot/differential/differential_abundance.py:265,276-277), and pooled uses the same expression through _auto_ls (kompot/differential/differential_expression.py:30,55-57). Measured 41.72104835510254 for both. On the constructor side, sync_parameters=True hands the same d, mu and ls to both DensityEstimators, while False hands neither any of them, so each estimates its own. That is exactly "separate".

• da() silently ignores two GPSettings fields. kompot.da() reads only n_landmarks, landmarks, ls_factor, batch_size, jit_compile and random_state from gp (kompot/anndata/differential_abundance.py:104-109). Measured against the default run, max|lfc − default| is:

  • GPSettings(ls=0.05): 1.8e-15
  • GPSettings(ls_scheme="pooled"): 1.8e-15
  • the same ls=0.05 passed through **density_kwargs: 1.65 (the positive control)

  sync_parameters itself is not named anywhere on da(), in its docs or in the DA CLI keys (kompot/cli/da.py:295-301). It reaches the model only by falling through **density_kwargs (line 202), where it does take effect: 3.0e-2. So the thing to integrate is da() honouring GPSettings, not a new parameter. The implementation should make da() honour or reject every GPSettings field, not only the new one.

• Not every value carries over to da yet. separate and pooled are exact translations of today's behaviour. condition1 and condition2 generalise directly, since each just estimates d, mu and ls from one condition. symmetric is a geometric mean of two length scales (differential_expression.py:73) and has no definition for d or mu. I would have da reject it with an explicit error until it has a definition and a measurement, rather than invent one.

• Swap behaviour differs by source, and the docs should say so per value. symmetric is swap-invariant by construction: the shared ls is bit-identical when the inputs are exchanged (49.74640341937494 in both orientations). pooled is not, because it runs an approximate nearest-neighbour search over the stacked union, whose row order depends on orientation: 41.7210 vs 41.7112, a relative gap of 2.4e-4 at this size. compute_d_factal subsamples 500 query cells above 500 combined cells, so under pooled d is order-dependent as well: 1.98779 vs 1.98875 at 1,000 cells.

• Out of scope for the scheme, and the same shape of defect: automatic landmarks come from the stacked union regardless of scheme, as discussed above. There is also the entry-point default split, where da() uses n_landmarks=None without gp and 5000 with one. It is the same problem as the ls_scheme default: a GPSettings default that is really a per-entry-point default. The None sentinel fixes both the same way, so I'd do them together.

Scope: this is a design answer, not a change to this PR, which remains documentation corrections on the length-scale feature. If this shape works for you, I can open a dedicated issue for the implementation: the GPSettings field, da() honouring gp.ls and the scheme, the sync_parameters deprecation alias, and the CLI keys on both sides. The one decision to make before the next release is the name, since ls_scheme ships under whatever name it has then.

[settylab-dotto-bot]

Skeptic verdict on f72573e, the corrected body, and the retraction: CHECK

The corrections are right, and every number in them re-measures. There is one scope gap in the new remedy sentence, and three presentation defects. None is a wrong measurement.

Verified (my own arms, run with random_state at its default. Metric: max relative gap fit(X1,X2) vs fit(X2,X1), with signed quantities negated):

arm	combined cells	log_fold_change	density
control: class defaults	400+400	0	0
sync_parameters=True, no landmarks (reference)	400+400	5.1e-01	2.767e-03
same, one shared landmark array	400+400	5.5e-01	2.771e-03
mu+ls explicit	400+400	5.6e-01	2.834e-03
d+mu+ls explicit	400+400	0	0
mu+ls explicit	200+200 / 250+250 / 150+350	0 / 0 / 0	0
mu+ls explicit	251+250 / 150+650	2.4e-01 / 1.9e+00	1.5e-04 / 3.0e-03
d+mu+ls explicit	150+650	0	0

(log_fold_change max-relative is inflated here by cells where it crosses zero. The zero/non-zero pattern is what the claims rest on.)

• (a) holds. Sharing landmarks leaves the sync half where it was (2.771e-03 vs 2.767e-03). The fit log shows why: ls itself differs between orientations (12.8327 vs 12.8268), not just d.
• (b) The 500 boundary is exact and correctly worded. mellon's compute_d_factal subsamples only if n < x.shape[0], where x is the stacked union. So the cutoff is strictly more than 500, counted over combined cells rather than per condition. 500 exactly is exact, 501 is not, and an unequal 150+350 is exact. "mu+ls suffice below 500" is now measured at three sizes, not just argued.
• The AST-identity claim holds. tests/test_da_swap_equivalence.py at c907da9 vs f72573e: the ASTs are equal once module, class and function docstrings are stripped. Positive control: the unstripped ASTs differ, so the comparison can see a difference.
• The retracted wording survives nowhere. Across the CHANGELOG, the PR body, the module docstring and the test docstring, the remedy is scoped to the landmark half and the corrected sentence is present at all four sites.

Findings:

1. "Passing all three is always correct" is true only for the sync_parameters half, and the configuration where it fails is the one the same paragraph calls the common one. With sync_parameters=True plus automatic landmarks, passing all three stays approximate. The measurement is d+mu+ls with n_landmarks=40 at 400+400: uncertainty gap 8.4e-01, density 5.2e-04. With one shared landmark array plus all three, the gap is 0. The two remedies compose, and neither suffices alone. This configuration is reachable from the documented entry point: da(adata, gp=GPSettings(...), sync_parameters=True) gets automatic landmarks, and sync_parameters reaches fit through **density_kwargs (established by reading the code, not measured through da()). A user following the sentence literally there is not exact. Fix: at all four sites (the retraction comment carries the same sentence), write "with sync_parameters=True, passing all three is always correct. If landmarks are automatic, also pass one landmarks array: both remedies are needed together."
2. The CHANGELOG is no longer concise, and the remedy exists only there. [Unreleased] went from 49 to 65 lines. The DA bullet is now 24 lines in three paragraphs, carries a raw probe number (4.784e-03 vs 4.720e-03) and editorial ("arguably the one you are least likely to be in"). None of the 0.8.0 entries does any of this; their long form is reserved for user-visible impact. Meanwhile no library docstring carries the precondition: the sync_parameters and landmarks entries in DifferentialAbundance.fit are silent on argument order, and sync_parameters even says it ensures "consistent density estimation across conditions". Fix: do what was already done for DE. Put the precondition and both remedies in DifferentialAbundance.fit's sync_parameters/landmarks docstrings, and cut the CHANGELOG bullet to about four lines that point there.
3. The PR body's suite section is stale and garbled. It still says the whole-suite run including the 8 DA tests is "outstanding". That run happened: 2131 collected, 1 failed / 2074 passed / 56 skipped, with the DA file 8/8 and the same single pre-existing test_dm_missing_pca failure. Just before that, a sentence from an earlier revision survives spliced in ("…not merely the count. The +13 passed is exactly the 13 … and the failing-ID set is identical, not merely the count"). Fix: replace it with the final totals.
4. The body's two numbers for sync_parameters=True without landmarks come from different datasets, and neither is labelled. The "Two configurations are only approximate" table (1.220e-02) sits directly under a paragraph naming 400+400 in 2-D. The archived probe behind it actually used 120+40 cells in 3-D, which is below the 500-cell boundary this same section then introduces. The 4.720e-03 "reference" a few lines later is presumably the 400+400 run. The same configuration therefore shows two different numbers with no stated reason. Fix: name the dataset beside the table.

Scope note: only the corrections were reviewed. The API-design question is being answered separately.

[settylab-dotto-bot]

Correction to the retraction above — c962036

The retraction comment I posted earlier says "Passing all three is always correct." That is
wrong as written, and a review pass on the corrections themselves caught it. I am not editing that
comment, because an edited comment hides that the claim was ever made.

It is true of the parameter half only. Passing d, mu and ls explicitly fixes the
parameters and leaves the landmark half untouched. Measured at 400 + 400:

configuration	uncertainty, between orientations
sync_parameters=True + automatic landmarks + d,mu,ls explicit	8.4e-01
sync_parameters=True + shared landmarks + d,mu,ls explicit	0

The two remedies compose and neither substitutes for the other. With sync_parameters=True,
passing all three is always correct — and if landmarks are automatic, pass one landmarks array
as well.

Corrected in the PR body, CHANGELOG.md, and tests/test_da_swap_equivalence.py. The precondition
and both remedies now also live in DifferentialAbundance.fit, under sync_parameters and
landmarks, which said nothing about order dependence before — that is where a user meets the
problem, and the CHANGELOG bullet is now a nine-line pointer rather than a transplanted essay.

---

This is the fourth instance on this PR of one habit, and the first where the over-generalised
measurement was not my own: I transcribed a number taken with landmarks controlled and published
it unconditionally. Re-reading it would not have helped — the condition a measurement was taken
under is not visible in the number, which is why each of these has been caught by someone
re-measuring rather than re-reading, and never by the person who made the claim.

The other three, for the record: "at its defaults" (true of the class, not the entry point users
configure through), "shared landmarks restore exactness" (true of the landmark half, not the
sync half), and "d is identical" (true at ≤ 500 combined cells, not above).

Also corrected in the same commit: the suite section quoted a superseded run and called the
abundance whole-suite check outstanding — it has since run, 1 failed / 2074 passed / 56
skipped, the single failure pre-existing and identical to base — and the 1.220e-02 comparison
table was unlabelled and read as the 400 + 400 run when it is 120 + 40 in 3-D.

Tracking: settylab/dotto-nexus#310.

[settylab-dotto-bot]

Skeptic re-check of fixes 1–4 at c962036: CREDIBLE

All four fixes are present and correct, and the new composed remedy survives the configurations it had not been measured on.

• F1 is fixed at every site. The CHANGELOG, the PR body, both test docstrings and the follow-up comment all say with sync_parameters=True, pass d/mu/ls, and if landmarks are automatic also pass one landmarks array. Scope attack on the composed sentence. Sync, shared landmarks and explicit d/mu/ls together mirror exactly (0 on log_fold_change, uncertainty and density) in four new settings: unequal 150+650; 3-D at 1500+1500 with 100 landmarks; landmarks = all 800 points (the DA: uncertainty is numerically unstable when n_landmarks >= n_combined (the GPSettings default on <5000 cells) #32 regime); and with sample indices (the sample-variance path). Each arm has a control that reads non-zero: the same data without explicit parameters, or with d omitted, gives gaps of 2.7e-03 to 7.7e-03. So the sentence needs no size or landmark-count qualifier, and structurally it cannot: at c962036, fit touches the stacked union in exactly two places, landmark selection and the sync block.
• F2 is fixed. The precondition and both remedies are now in DifferentialAbundance.fit under landmarks and sync_parameters, and the CHANGELOG bullet points there. [Unreleased] is 50 lines. Nit, not blocking: the landmarks entry says "with random_state=None the selection is not even reproducible". None is the class default, so saying "(the default)" would make that land.
• F3 is fixed. The suite section reads 1 / 2074 / 56 at f72573e, the +26 is accounted for, and the spliced sentence is gone. Not re-running the suite is sound: differential_abundance.py and test_da_swap_equivalence.py are AST-identical to f72573e with docstrings stripped (the raw ASTs differ, as the control), and pytest collects no doctests (doctests = true in pyproject.toml belongs to flake8).
• F4 is fixed. The table is labelled "120 + 40 cells in 3-D".

[katosh]

@settylab-dotto-bot fully implement the suggestion with a new worker and make sure everything is complete, consistent and we resolve all remnants of the previous implementation.

[settylab-dotto-bot]

Skeptic findings fixed: 298e032

Must-fix, mu at a pinned d. With d pinned alone, "condition1", "condition2" and "symmetric" shared a mu computed at the fractal d. A pinned d now feeds mu under every scheme. For "symmetric", it also replaces the pooled d. "pooled" was already correct.

New tests pin d alone and ls alone for each of the three schemes. Each test's reference is the estimator's own derivation, not the scheme under test. Run against e9f798b's library code, the three d-alone tests fail and pass at 298e032. The three ls-alone tests pass on both, as controls, since ls does not enter d or mu.

Minor (a). Passing param_scheme= as a flat keyword to da() or de() now raises a ValueError that points at gp=GPSettings(param_scheme=...). Previously it raised a TypeError about multiple values for one keyword. The new test fails at e9f798b and passes here.

Minor (b). The CHANGELOG argument-order bullet moved from Added to Documentation.

Tests. The four touched test modules and all seven test_anndata* suites give 211 passed and 1 skipped at 298e032. I did not rerun the full suite. The one change outside differential abundance and settings is a guard in de() that fires only when param_scheme is passed as a keyword.

[settylab-dotto-bot]

Skeptic re-check: CREDIBLE at 298e032, the remote tip, re-checked before posting. The re-check was scoped to the three items from my previous verdict, and all three are resolved.

1. A pinned d now feeds mu. I re-ran my original probe on the same data (200 + 160 cells, 4-D, d=2.0):
   • condition1 shares (2.0, -12.2865), which matches "separate" for condition 1 at the same pin (the estimator's own answer). At e9f798b this was -12.6685.
   • condition2 matches condition 2's own estimate at the pin.
   • symmetric shares d = 2.0, and its mu equals compute_mu over the pooled within-condition distances at d = 2.0.
   • The mirror and symmetric's swap bit-identity still hold with d pinned.
   • The new tests, run by me against both library trees (git archive, the 298e032 test file copied in): the three …derives_mu_at_a_pinned_d tests and the flat-keyword test fail on e9f798b and pass on 298e032. The ls-alone controls pass on both, as they should, which confirms that the partial-pin tests discriminate.
2. Flat param_scheme= is refused clearly. da(..., param_scheme=…), de(..., param_scheme=…) and the deprecated compute_differential_abundance(..., param_scheme=…) all raise a ValueError naming gp=GPSettings(param_scheme=...). compute_differential_abundance(sync_parameters=True) still works.
3. CHANGELOG. The argument-order bullet is now under Documentation.

Regression: test_param_scheme.py together with test_da_swap_equivalence.py at 298e032 gives 81 passed. I did not re-run the full suite: outside da, the delta is a four-line guard in de().

[katosh]

@settylab-dotto-bot we don't need to state "DE also" as this is the expected case and only the opposite would need mentioning. Why is param_scheme: str = "separate", the default in one of the da helpers? Make sure this PR does not change the default behavior of kompot . Can we drop sync_parameters entirely? No backwards compatibility or deprication warning needed at this point.

[settylab-dotto-bot]

All four points are addressed in 7e37ac6..812c15b.

Why param_scheme: str = "separate". It was in _check_da_overwrites, the helper that decides whether a da() run may overwrite an earlier one. I gave it a default so that the six direct calls in tests/test_da_helpers_overwrite.py kept working without the new argument. da() always passes the resolved value, so the default never ran, but it was a second copy of DA's default that would have drifted silently if the default ever changed. The DE helper _check_overwrites had the same pattern with "condition1". Both now take param_scheme as a required argument. A run recorded before param_scheme existed now reads back as DEFAULT_PARAM_SCHEME, so each entry point's default is defined in one place.

Default behaviour is unchanged. I compared 19ee1a1, the PR's merge base on main, with 812c15b. Each tree was exported with git archive, selected by PYTHONPATH (the run asserts kompot.__file__), and both ran in one Slurm job. Every entry point was called with no arguments beyond the data. Data: synthetic, 260 + 190 cells, 4-D, 25 genes, seed 0. np.random.seed(0) was set before each call so the automatic landmarks are reproducible without passing random_state.

entry point, at defaults	arrays	bit-identical	max abs diff
de()	9	9	0
da()	6	6	0
da(gp=GPSettings())	6	6	0
DifferentialExpression().fit + predict	7	7	0
DifferentialAbundance().fit + predict	7	7	0
compute_differential_expression()	9	9	0
compute_differential_abundance()	6	6	0
CLI kompot de, flags only	10	10	0
CLI kompot da, flags only	7	7	0
CLI kompot da -c with its shipped template	7	7	0
CLI kompot de -c with its shipped template	outcome	identical	fails on both trees, see below
total	75	75	0

Controls:

• Determinism: a second base run in the same job matches the first on 75 of 75 arrays.
• Positive controls differ from the base default: param_scheme="symmetric" in de() differs on 8 of 9 arrays (max 0.12), "pooled" in da() on 5 of 6 (max 0.27), and a relative 1e-7 change to ls_factor on 8 of 9 (max 2.1e-07).

The only difference at the defaults is metadata: the run record in uns gains gp.param_scheme ("condition1" for de(), "separate" for da()).

sync_parameters is removed entirely. There is no alias, warning or shim. param_scheme takes its argument position in DifferentialAbundance.fit, so allow_single_condition_variance stays where it is on main.

$ git grep -n sync_parameters 812c15b
812c15b:CHANGELOG.md:16: - `sync_parameters` on `DifferentialAbundance.fit` and `da()`; use `param_scheme="pooled"`.

A caller who still passes it gets Python's ordinary error, measured at head through DifferentialAbundance.fit, da() and compute_differential_abundance(): TypeError: DensityEstimator.__init__() got an unexpected keyword argument 'sync_parameters'. There is one exception, which is not specific to this keyword: when both density predictors are injected, fit builds no estimator, so it never checks any density keyword. This also settles the earlier warnings vs logger question, because this PR no longer contains a deprecation.

Wording. Docstrings, docs/source/cli.rst and the CHANGELOG now describe the expected case once and name only where DA differs. For example, the CLI table lists all GPSettings keys and adds "In DA, sigma is the density estimator's and eps and use_empirical_variance do not apply." The CHANGELOG [Unreleased] block is down from 43 lines to 25.

Tests. Full suite, both from branch-checked-out worktrees: 298e032 gives 2 failed, 2108 passed, 56 skipped; 812c15b gives 2 failed, 2105 passed, 56 skipped. The failing-ID sets are identical (comm -3 is empty): the pre-existing test_cli.py::TestCLIDiffusionMaps::test_dm_missing_pca and the worktree-only test_provenance.py::test_source_checkout_resolves_sha_and_editable (#33). The 3 fewer passes are the removed sync_parameters tests.

Pre-existing and out of scope: at base and head alike, kompot de -c de_config_template.yaml fails with DifferentialExpression.fit() got multiple values for keyword argument 'landmarks'. The DE CLI does not route the template's landmarks: null into GPSettings.

[settylab-dotto-bot]

kompot-pscheme-followup

I pushed the operator's four follow-ups on kompot PR #22 as 298e032..812c15b.

Full report: https://github.com/settylab/dotto-nexus/blob/d7f29608ecdcb360c2f845a8a36a541e09f55dc3/assets/22/pscheme2_2026-09-23_112528_kompot-pscheme-followup.md

[settylab-dotto-bot]

#32 is fixed on this PR at 6606ee2, and the fix knowingly changes da() results in one regime. My earlier "no default behaviour change" was stated for 812c15b and now holds only outside that regime; the scope is below.

Cause

• mellon, belongs upstream. Without ADVI, DensityEstimator's uncertainty is a diagonal Laplace approximation in whitened coordinates z (f = L z, L = chol(K_uu)). A diagonal covariance in z is not invariant to the order of the landmark rows, and at kompot's ls_factor=10, K_uu is near-singular, which amplifies the effect. In a mellon-only test the landmark set is the same and only its row order changes. On 400 cells:

  ls_factor	cond(K_uu)	uncertainty, rel. diff.	mean, max rel. diff.
  1	4.7e10	median 9e-6	2e-11
  10	1.4e16	median 0.94	5e-12

  This also means the part of DA: uncertainty is numerically unstable when n_landmarks >= n_combined (the GPSettings default on <5000 cells) #32 about "K_uu ill-conditioned in the sparse path" is not the mechanism. The conditional-variance term is about 1e-7. Almost all of the uncertainty comes from the Laplace term.

• kompot, fixed here. At n_landmarks >= n_combined, DifferentialAbundance.fit used the stacked union [X_condition1; X_condition2] as the landmark array. Its row order depends on which condition is passed first. Each condition's model therefore received the same landmark set in a different order, and the uncertainty differed.

Fix

When n_landmarks >= n_combined, no landmarks are built, and each condition gets its full GP on its own cells, the same as n_landmarks=None. No approximation was requested, and the cost is lower: n_c³ per condition instead of n_combined³.

What changes, entry point by entry point

entry point, at its defaults	in the #32 regime?	result
da() with no gp	no (n_landmarks=None)	bit-identical
da(gp=GPSettings(...)), ≤ 5000 cells	yes	changes, see below
da(gp=GPSettings(...)), > 5000 cells	no (k-means)	unchanged
DifferentialAbundance()	no (n_landmarks=None)	bit-identical
compute_differential_abundance()	no	bit-identical
CLI kompot da, flags only or the shipped template	no	bit-identical
CLI kompot da with any GP key but no n_landmarks (e.g. --random-state 0), ≤ 5000 cells	yes	changes
everything in de()	not affected	bit-identical

Inside the regime, measured on the #32 fixture (400 + 400 cells, da(gp=GPSettings(random_state=0)), 812c15b against 6606ee2):

• z-scores and PTPs move to the full-GP values. z changes by a median of 34% (q90 81%).
• Orientation disagreement falls to exactly zero. In z it was a median of 82% and is now 0. In PTP it was a median of 96% and is now 0. Fold changes already agreed and still do.
• Point estimates also move, slightly. Fold changes and log densities change by at most 6.0e-4 absolute (median relative 2e-5). The old model placed latent values at the other condition's cells too, so the result is close to, but not exactly, the old one.
• The new result matches the gp=None full-GP run to 3e-14, and is bit-identical to it at class level.

Outside the regime, 19ee1a1 against 6606ee2 in one job: 70 of 75 arrays are bit-identical. The 5 that differ are all from da(gp=GPSettings()) on 450 cells, which is in the regime: z up to 6.3e-2, PTP up to 3.7e-2, fold changes and densities up to 4.7e-6. Direction labels are unchanged. A control at n_landmarks=100 (below 800 cells) is bit-identical between 812c15b and 6606ee2 single-threaded, 10 of 10 arrays. Multi-threaded, the k-means landmark path varies by up to about 3e-10 from one process to the next on either tree; that variation predates this PR.

The CLI row with a GP key is measured: with --random-state 0, the orientation gap in z goes from 0.12 to 0.

Tests

test_landmarks_covering_every_cell_fit_the_full_gp[above|equal] checks two things at n_landmarks=5000 and at exactly 800: the output is bit-identical to n_landmarks=None on six outputs, and it is orientation-symmetric. It fails at 812c15b and passes here. A control confirms that 799 landmarks still builds 799.

Mutants, with predictions written first. All four came out as predicted:

• >= changed to > breaks only [equal].
• Counting only condition 1's cells breaks only the control.
• Dropping the explicit n_landmarks=0 passed to mellon survives. It only matters for a single condition above 5000 cells, and I did not add a 5000-cell full-GP test. That is a stated coverage limit.

Full suite at 6606ee2: 2 failed, 2108 passed, 56 skipped. The failing-ID set is identical to 298e032: test_dm_missing_pca and the worktree-only provenance test (#33). There are 3 more passes than at 812c15b, which are the new tests. 794bf7f changes only .github/workflows/tests.yml.

Not fixed here

With fewer landmarks than cells, the uncertainty still depends on argument order, and by as much. On this fixture, with k-means at 100 landmarks, it differs by a median of 87%. This is the same mellon mechanism, applied to k-means centroids whose order also depends on the input. Passing one landmarks array to both runs makes it exact, which is already documented and tested. The real fix, an uncertainty that does not depend on landmark order, belongs in mellon. I have drafted a mellon issue with a mellon-only reproduction but not filed it; it needs your go-ahead.

CI, separate commit 794bf7f: this bounds and diagnoses the intermittent 3.11 hang. It sets a 60-min job timeout and fail-fast: false, and adds pytest-timeout at 300 s with the thread method, which dumps each thread's stack when a test hangs, plus --durations=25. It limits the hang and captures evidence; it does not cure it, and the root cause is still being hunted.

[settylab-dotto-bot]

Skeptic verdict: CREDIBLE at 794bf7f (reviewed against main = 19ee1a1, which is also the merge base). The results change only where the CHANGELOG says, which is da(gp=GPSettings(...)) when n_landmarks is at least the number of cells. sync_parameters is gone and fails loudly, and there is now one default per entry point. There is one small API break to fix before merge (finding 1).

I used my own harness: git archive of both trees, each selected by PYTHONPATH with kompot.__file__ asserted, all in one Slurm job. The inputs are synthetic and my choice, not defaults: 280 + 210 cells, 5-D, 30 genes, and np.random.seed(0) before every call.

Default behaviour

• Small data (490 cells): 207 of 212 arrays are bit-identical between main and 794bf7f, and main against itself is 212 of 212. I ran 18 setups:
  • de(), da(), both with GPSettings(), both with sample_col, and both with 60 landmarks;
  • the compute_differential_* wrappers, and both classes via fit/predict;
  • a second run on the same object, and a rerun on an h5ad saved by main;
  • the CLI with flags only and with each tree's template.
• The 5 arrays that differ are all da(gp=GPSettings()), which is in the DA: uncertainty is numerically unstable when n_landmarks >= n_combined (the GPSettings default on <5000 cells) #32 regime. There, fold changes and densities move by at most 5e-6, z by at most 0.079 and PTP by at most 0.093. No direction label changes.
• Controls differ as they must. de with symmetric differs on 8 of 13 arrays, da with pooled on 5 of 10, and a 1e-7 ls_factor nudge on 13 of 13.
• The saved run record gains only gp.param_scheme, which stores the resolved value. Logs and warnings are identical apart from timestamps.
• Large data (5,700 cells, true defaults, outside the regime): 23 of 23 arrays are bit-identical, run single-threaded on one node, and main against itself is 23 of 23. Multi-threaded across nodes, main differs from itself at about 2e-12 in DE, so only a single-threaded run gives an exact comparison.
• Which entry points are in the regime. da() without gp and DifferentialAbundance() are unchanged, matching your correction. The DE template's landmarks collision fails identically on both trees.

The #32 fix

• On da(gp=GPSettings(random_state=0)) with 490 cells, I compared forward z against minus the reverse z. The gap is 0.214 (max) on main and exactly 0 at head. For log fold change it is 4.6e-14 on main and exactly 0 at head.
• At head, da(gp=GPSettings()) equals da() to within 2.4e-14. The only difference between the two paths is batch_size (100 against None), so this is roundoff, and the fix does produce the full GP.
• Your mutation record matches its predictions, which were written first. M3 surviving is explained correctly.

sync_parameters

• git grep finds only the CHANGELOG line.
• It now raises an error through DifferentialAbundance.fit (by keyword, True or False), through da(), through compute_differential_abundance() and through a CLI config key. Passed by position, True or False, it raises a ValueError. All of these ran without error on main.
• Not a blocker, and not new: with density_predictor1/2 injected, fit builds no estimator and silently accepts any keyword. I checked this with totally_bogus=1 on both trees.

Other checks

• Helper default: no hard-coded scheme default remains. Each module has one DEFAULT_PARAM_SCHEME, and both overwrite helpers import it.
• Wording: no "DE also"-style cross-references remain in the diff, the docs or this PR's body.

Findings

1. Fix before merge, a one-line change. DifferentialExpression.fit adds param_scheme as the 9th parameter, ahead of landmarks. So fit(X1, Y1, X2, Y2, 1.0, None, 10.0, L) works on main and fails at head with "The truth value of an array with more than one element is ambiguous". Moving param_scheme after allow_single_condition_variance keeps positional calls working. In DifferentialAbundance.fit it took sync_parameters' position, which is fine.
2. Docs, minor. The n_landmarks docstring and the CHANGELOG say the regime now fits the full GP "as with n_landmarks=None". That holds only up to 5,000 cells per condition. Above that, None leaves mellon at its default of 5,000 landmarks per condition, while the explicit n_landmarks=0 fits the full GP. This is not a regression, because the old code used the whole union as landmarks there. Suggested wording: "each condition gets its full GP".

Regression tests. I ran the full suite on detached worktrees and compared the sets of failing test IDs:

• 19ee1a1: 1 failed, 2028 passed.
• 812c15b: 1 failed, 2106 passed.
• 794bf7f: 1 failed, 2109 passed. The 3 new DA: uncertainty is numerically unstable when n_landmarks >= n_combined (the GPSettings default on <5000 cells) #32 tests pass.

On each tree the only failure is test_cli.py::TestCLIDiffusionMaps::test_dm_missing_pca, which already fails on main.

[settylab-dotto-bot]

kompot-pscheme2-skeptic

Skeptic verdict credible at 794bf7f on kompot PR #22, reviewed against main 19ee1a1. - Default behaviour is unchanged except inside the #32 regime.

Full report: https://github.com/settylab/dotto-nexus/blob/56c84186cf60a175359f73e5bdf70e84030f6bbf/assets/22/pscheme2sk_2026-09-23_115803_kompot-pscheme2-skeptic.md

[katosh]

@settylab-dotto-bot fix the last failing test, file the Mellon issue for our consideration, investigate implications on accuracy and runtime through our benchmarks.

[settylab-dotto-bot]

Follow-ups: the skeptic's two items (31315ea, 343486f) and your comment (da733b4, settylab/Mellon#22).

• Positional order restored. param_scheme is now the last parameter of DifferentialExpression.fit, so fit(X1, Y1, X2, Y2, 1.0, None, 10.0, L) works again. Checking all 56 public callables against main found one more shift: GPSettings, where GPSettings(1.0, None, 10.0, 100) now binds n_landmarks again. The only remaining difference is intended: DifferentialAbundance.fit puts param_scheme where the removed sync_parameters was, so allow_single_condition_variance keeps its position, and a stale positional bool raises. Two new tests call both positionally in main's order. Both fail at 794bf7f and pass here. The fit test fails with exactly the reported error; the GPSettings test fails with param_scheme ... (got 100).
• Docs. The full-GP fallback equals n_landmarks=None only up to 5000 cells per condition, and the CHANGELOG and docstrings now say so.
• test_dm_missing_pca is fixed in the code, not the test. kompot dm imported Palantir before validating its input, so without Palantir installed a missing PCA key was never reported. CI installs .[all], so it passed there. The input is now checked first. The test fails at 794bf7f and passes here, and a new test covers the Palantir-missing error when the input is valid.
• Mellon issue, for consideration: DensityEstimator uncertainty depends on the row order of landmarks (for consideration) Mellon#22. It has a minimal reproduction: the same landmark set in two row orders changes the uncertainty by a median of 94% at ls_factor=10, while the mean is unchanged to 5e-12. It also states what I measured and what I did not.

Full suite at da733b4: 1 failed, 2112 passed, 56 skipped. The only failure is the worktree-only provenance test (#33), which does not fail in CI's checkout. test_dm_missing_pca has left the failing set.

[settylab-dotto-bot]

Skeptic re-check: CREDIBLE at da733b4. All four items hold.

1. Positional order. I compared the signatures of all 192 public callables across every kompot submodule, 19ee1a1 against da733b4, using my own inspect pass.
   • DifferentialExpression.fit and GPSettings now only add param_scheme at the end, so every existing positional prefix is unchanged.
   • The only shift left is DifferentialAbundance.fit, where param_scheme replaces sync_parameters, as accepted.
   • Behaviour check: fit(X1, Y1, X2, Y2, 1.0, None, 10.0, L) uses L as landmarks on both main and da733b4. At 794bf7f it failed. GPSettings(1.0, None, 10.0, 123) sets n_landmarks=123 again.
   • No kompot code constructs GPSettings positionally.
2. Docs. The CHANGELOG, the n_landmarks docstring and the code comment now say "equals None only up to 5000 cells per condition". That matches mellon's DEFAULT_N_LANDMARKS = 5000 and its full-GP rule.
3. test_dm_missing_pca. The fix only moves the Palantir import below the input checks; there is no numerical change. I ran da733b4's tests from a neutral directory, so PYTHONPATH selected the library, without Palantir installed:
   • 794bf7f library: test_dm_missing_pca, test_de_fit_keeps_the_positional_order_of_main and test_gpsettings_keeps_the_positional_order_of_main fail.
   • da733b4 library: 4 passed, 2 skipped (they need Palantir).
4. DensityEstimator uncertainty depends on the row order of landmarks (for consideration) Mellon#22. I ran the snippet exactly as posted. It prints both output lines exactly as they appear in the issue. The issue contains no internal identifiers.

[settylab-dotto-bot]

Benchmark check of this PR: accuracy and runtime, main @ 19ee1a1 vs PR head @ 794bf7f

The head has since moved to da733b4. Relative to 794bf7f that adds docstrings, moves param_scheme within the DifferentialExpression.fit / GPSettings parameter order (keyword calls are unaffected, and every call here uses keywords), and adds kompot dm CLI input validation. None of it is reachable from these runs, so the results below carry over unchanged.

Both kompot trees were installed from git archive exports into one fixed environment (Python 3.13.5, jax 0.9.2, numpy 2.4.4, mellon 1.7.1). All arms of each comparison ran back to back on the same Slurm node.

1. Reference benchmarks: no output changes

Every benchmark kompot call goes through de() with param_scheme=None, which resolves to "condition1": the old code path, with no extra computation and the same length-scale inheritance. The only da() calls pass n_landmarks=None and no ls, so neither DA change in this PR can reach them. There is no DA accuracy benchmark. Measured:

• Spike-in DE benchmark (7 scenarios, 16,285 genes, fixed seeds): the kompot score is bitwise identical main vs PR in every scenario. So are all 7 per-gene kompot outputs and every other method in the run. This holds at two spike-set sizes (150 and 500 positives). A main-vs-main repeat is also bitwise identical, so the noise floor is exactly 0. AUROC/AUPRC are therefore unchanged (e.g. 0.9612 / 0.7157 on the uniform scenario, both trees). One scenario would not build at the 500-positive setting and is covered by the 150-positive runs.

• Runtime benchmark (kompot DE, wall time, 3 repeats per tree, interleaved):

  cells × genes	main (s)	PR (s)	Δ	within-tree repeat spread
  5,000 × 5,000	64.1	64.2	+0.1	4.1 / 4.4
  20,000 × 5,000	112.9	113.2	+0.3	1.5 / 8.2

  Peak process memory is identical (17.4 / 21.3 GB).

2. Targeted probe of the #32 regime (not a benchmark)

The benchmarks don't cover the path this PR changes, so I probed it directly: da(gp=GPSettings(random_state=0)) on a real 40-dim diffusion embedding, balanced subsamples of two conditions, both orientations, 2 repeats. Plain da() (exact full GP) serves as the reference.

cells (combined)	max |Δz| main vs PR (fwd / swap)	max |Δlfc|	orientation gap max|z_fwd + z_swap| main → PR	wall s main → PR	peak RSS main → PR
500	0.21 / 0.38	2.5e-5	0.54 → 0	31.2 → 29.6	1.24 → 1.17 GB
2,000	1.41 / 8.00	1.3e-5	9.41 → 0	36.0 → 32.5	1.96 → 1.46 GB
4,800	2.51 / 9.79	1.9e-3	10.89 → 0	64.5 → 43.1	2.77 → 1.57 GB
5,800 (control, above 5,000)	≤ 2.7e-6	≤ 5.6e-7	6.10 → 6.10	74.0 → 73.1	3.08 → 3.05 GB

• Below 5,000 cells, the PR's da(gp=GPSettings()) equals plain da() to ≤ 9e-14 in z. The z shifts in the table are the error on main, which this PR removes. The fraction of cells with |z| > 2 at 4,800 cells goes from 0.341 / 0.320 (fwd / swap) to 0.335 in both orientations.
• Noise floor: results are deterministic across repeats, except on the >5,000 landmark path, where |Δz| ≤ 2.7e-6. The 5,800-cell control's main-vs-PR delta sits at that floor, so it is unchanged. Plain da() is bitwise identical main vs PR at every size. Plain-da() wall times agree between trees within 0.7 s. Wall times include roughly 25 s of JIT compilation.
• Point estimates move more at larger sizes than the PR body's measured cases (≤ 4.7e-6 on 450 cells, ≤ 6.0e-4 on 400 + 400): here |Δlfc| is 2.5e-5 at 500, 1.3e-5 at 2,000 and 1.9e-3 at 4,800 cells (2,400 + 2,400). That is still small next to the z shifts.
• The orientation gap above 5,000 cells (6.10 on both trees) is the pre-existing landmark order dependence the PR already documents. It is unchanged here.

Independent validation of these numbers is pending. Corrections, if any, will follow in this thread.

[katosh]

@settylab-dotto-bot runtime reports without JIT please, and assess prediction accuracy on the semisynthetic spike in data.

[settylab-dotto-bot]

Follow-up: runtime without compilation, and accuracy against spike-in ground truth (main @ 19ee1a1 vs PR @ 794bf7f; da733b4 adds only docs, parameter order and CLI validation)

Correction to my previous comment: "one scenario would not build at the 500-positive setting" was wrong. I used a retired setting for that scenario. With the correct setting it builds and is bitwise identical main vs PR (AUROC 0.9762 / AUPRC 0.6658 on both). So all 7 spike-in scenarios are bitwise identical at both spike-set sizes. "Bitwise" means same node throughout: see the noise floor note at the end.

1. Runtime without compilation

GPSettings.jit_compile already defaults to False, and da() never forwards it to mellon, whose jit also defaults to False. Every earlier number was therefore already "jit off". The remaining ~26–29 s per process is JAX compiling each operation the first time it runs, which no flag disables. Measured two ways: cold = the first call in a fresh process, warm = an identical second call on the same shapes, so compilation is excluded. Exclusive node, trees interleaved, 6 runs per cell, times in seconds:

cells	cold main → PR	warm main → PR	compilation (cold − warm) main / PR
500	29.8 → 28.1	0.78 → 0.71	29.1 / 27.4
2,000	32.1 → 29.7	2.76 → 1.65	29.3 / 28.0
4,800	48.6 → 36.0	20.4 → 8.5	28.2 / 27.4
5,800 (control)	55.6 → 55.7	26.7 → 26.8	28.9 / 28.9

Plain da() (same code on both trees) agrees within 0.13 s cold and 0.26 s warm, which gives the timing noise floor. The PR's speedup is in the computation, not the compilation: at 4,800 cells, warm time drops by 58%, and only 0.8 s of the 12.6 s cold saving is compilation. The earlier 64.5 → 43.1 s ran on a different CPU type and wasn't split into cold and warm. By this split, essentially none of that saving was compilation either.

The runtime benchmark (DE, which this PR doesn't touch) behaves the same: at 5k cells, 59.1 → 59.0 s cold and 39.5 → 39.3 s warm (compilation ≈ 33%). At 20k cells, 101.8 → 101.8 s cold and 97.3 → 95.0 s warm, where the 2.3 s gap is inside the 3.1 s within-tree repeat range.

2. Accuracy against spike-in ground truth

• What truth exists: the spike-in data only carries per-gene expression truth. Both conditions contain identical cells, so there is no abundance truth, and at full size the DA: uncertainty is numerically unstable when n_landmarks >= n_combined (the GPSettings default on <5000 cells) #32 path doesn't fire. The DE result is bitwise identical (above), so DE accuracy is unchanged.
• DA truth has to be planted. I used the same cells and embedding: subsample (1,200 or all 2,917 cells), split at random into two equal halves, then remove 50% or 75% of one cell-type region (two regions, ~9% of cells each) from one half. There is also a null arm with no removal, and 3 random seeds. Truth = in-region cells. Everything here is in the DA: uncertainty is numerically unstable when n_landmarks >= n_combined (the GPSettings default on <5000 cells) #32 regime. These are planted, illustrative shifts, not a sample of real biology.
• The PR's da(gp=GPSettings()) equals exact-GP da() to ≤ 5e-12 in z in every configuration. On main, accuracy depends on which condition is passed first:

cells	region / kept	AUROC main fwd	main swap	PR (both)
2,917	A / 25%	0.631	0.760	0.716
2,917	A / 50%	0.506	0.610	0.575
2,917	B / 50%	0.976	0.969	0.976
1,200	B / 25%	0.931	0.947	0.939

• Across 48 paired draws, main − PR ranges from −0.18 to +0.08 in AUROC and −0.39 to +0.22 in AUPRC (mean |Δ| 0.018 and 0.033). Averaged over both argument orders, main is within ~0.02 AUROC of the PR in every configuration. So the fix doesn't raise expected accuracy; it removes an order-dependent error of up to ~0.18 AUROC in a single run. The SD across seeds is 0.004–0.30 AUROC, so the main-vs-PR differences are paired within a draw.
• Calibration of z (pre-existing, unchanged by this PR): in the null arm (random halves, nothing planted), z carries a seed-dependent global offset, identical on both trees and on exact-GP da(). Mean z reaches −2.17 at 1,200 cells (66% of cells with |z| > 2) and −1.05 at 2,917 cells (21%); other seeds are close to 0. This fits the DA default param_scheme="separate" (per-condition d/mu/ls), but I haven't tested whether "symmetric" removes it. Power at |z| > 2 is low for the weaker region (≤ 0.19) and ~0.96 for the stronger one at 2,917 cells.

Noise floor, two levels: repeats on one node are deterministic (exactly 0 below 5,000 cells, ~1e-6 on the landmark path above it). Across node types, identical inputs give z values that differ by up to 0.2–0.8 on both trees, including exact-GP da(). Every main-vs-PR comparison above is paired on one node, and the main-vs-PR deltas replicate across the two node types (4,800 cells: Δz 2.51 / 9.79 vs 2.50 / 10.18; the PR's orientation gap is 0 on both).

[settylab-dotto-bot]

Independent validation: CREDIBLE. Both benchmark comments hold. I re-measured them (not re-read) on the current head da733b4 against main @ 19ee1a1. Four corrections are wording or attribution only, and none changes a conclusion.

Re-measured, in my own clone with my own drivers and comparison scripts, main vs da733b4 paired on one node:

• Reachability. I re-derived every kompot call site in the benchmark suite. All DE calls use keywords with no param_scheme. The two DA calls use the legacy wrapper with n_landmarks=None and no ls. No benchmark uses the CLI. Neither DA change can reach them.
• Spike-in DE benchmark, 500 positives, two scenarios (including the one the first comment wrongly reported as unbuildable): all kompot output columns are bitwise identical, as are the other methods in the run. My main arm also equals the earlier main arm on a different node.
  • Positive control: a planted 1e-6 relative change to one length scale altered 6 of 11 kompot columns (Mahalanobis max |Δ| 4.3e-5).
  • AUROC and AUPRC did not move at 6 decimals under that change. The bitwise check can see a change of that size and the metric cannot.
• Targeted DA: uncertainty is numerically unstable when n_landmarks >= n_combined (the GPSettings default on <5000 cells) #32 probe.
  • 4,800 cells: the PR equals exact-GP da() to 8.7e-14 in z, and main differs by |Δz| 2.51 / 9.79 (fwd / swap). The orientation gap is 10.9 on main and 0 on the PR.
  • 5,800-cell control: unchanged, with every tree delta ≤ 2.3e-6 in z, inside the repeat floor.
  • Plain da() and the legacy wrapper: bitwise identical.
• Speed and memory at 4,800 cells:
  • Cold time drops 53.0 → 37.5 s and warm 26.7 → 11.3 s (−58%). Compilation costs 26.3 s on both trees, so none of the saving is compilation.
  • Peak RSS drops 2.76 → 1.57 GB, and the 5,800 control is unchanged.
  • The mechanism is structural. On main, mellon uses all 4,800 combined cells as fixed landmarks for each 2,400-cell condition, where the PR fits the exact 2,400-cell GP.
• JIT, checked in source: da() stores GPSettings.jit_compile and never passes it to mellon, whose jit default is False. The earlier timings were already jit-off, as the follow-up says.

Accuracy (planted-truth design): fair, and correctly bounded. Main and PR get identical inputs, seeds and node. The spike-in data carries no abundance truth, because both conditions are the same cells, and the follow-up says so rather than overclaiming. From its per-draw records:

• No systematic bias. Averaged over both argument orders, main − PR is negative in 13 draws and positive in 11, with median −0.0001 AUROC. "Order dependence, not bias" is supported.
• The null-arm z offset is pre-existing. It is identical on the PR, on exact-GP da() and on main's exact GP.

Corrections, all minor:

1. Draw count. "48 paired draws" is 24 draws × 2 argument orders.
2. "Within ~0.02 of the PR in every configuration" holds for the seed-averaged rows. For a single draw, the order-averaged gap reaches 0.051 AUROC.
3. The single-order error is not random in direction. The depletion is always planted in the same condition. When that condition is passed second, main scores below the PR, in 6 of 6 draws in the weaker region at 2,917 cells. Passed first, it scores above.
4. The z noise floor follows the core count, not the node type. On one node, the identical exact-GP da() call run with 8 versus 4 visible cores changed z by up to 0.35. 8-core runs on two different nodes were bitwise identical. The follow-up's timing node exposed 36 cores.
   • This is a pre-existing numerical sensitivity of da()'s z, unrelated to this PR.
   • Every main-vs-PR comparison was paired within one job, so no conclusion changes.

The first comment's scenario error has been corrected explicitly in the follow-up. It is now verified: all 7 scenarios are bitwise identical at both spike-set sizes.