Sample variance: bound its memory, fix the planner, and document what it costs (0.9.0)

[settylab-dotto-bot]

Fixes #25
Fixes #26
Fixes #27
Fixes #23
Fixes #28
Fixes #30
Fixes #31
Fixes #33

Cuts 0.9.0. Three issues that are really one story: sample variance costs
much more than it needs to, the tool you would use to find that out was lying,
and nothing in the documentation said any of it.

It now also carries every other open issue: see "Bundled into 0.9.0" directly below.

Bundled into 0.9.0: every other open issue (added 2026-09-25)

The maintainer asked for every remaining open issue to ship in 0.9.0, so this
PR now also fixes #23, #28, #30, #31 and #33, and carries one
review change to the sample-variance guidance. main (#22, #32) is merged in;
the only conflict was CHANGELOG.md, and #22's [Unreleased] entries now sit
under 0.9.0. Each fix has a test that fails without it, verified against a
git archive of the pre-fix merge commit.

Revised after an adversarial review of a308f94 (533ab0f, b1a85f2). The
review found that the #30 check was blind to one leak shape and passed while 14
leaks were live, so the "clean" result below was the second false zero on
#30 in this PR. It also found that the #28 fix could stamp a wrong sha. Both
are fixed; details are in the #30 and #28 sections. Head after the revision: eac37df.

issue	disposition	commit	proof it bites
#23	fixed	74e62c6, 151f9b8	CLI suite on anndata 0.12.6: 49 failed → 0
#28	fixed (with #33)	04c8574, 533ab0f	6 new provenance tests fail on base; the 2 vendored-tree tests fail without the root guard
#30	fixed	da7979c, 533ab0f	build test fails with the hook neutered, and again with only the new pattern removed
#31	fixed (deprecated)	e6ddfc7	both behavioural tests fail on base
#33	fixed (same root as #28)	04c8574	as #28

#23: the CLI could not write under anndata 0.11 and 0.12

anndata 0.11 and 0.12 refuse to write a pandas nullable-string array unless
anndata.settings.allow_write_nullable_strings is set, and pandas 3 produces
one for every string column, the obs index included. anndata 0.13 defaults
the setting to None, which allows the write, and that is why CI (anndata
0.13.4) never saw it.

The 49 failures on the issue were not raised by the CLI. They came from the
tests' own fixtures writing their input with a plain write_h5ad, before any
Kompot code ran. The CLI would have failed next, at its output write, so both
halves were real:

• kompot.cli.utils.write_output is now the one writer for de, da,
  smooth and dm. It opts in with anndata.settings.override(...) for the
  duration of the write only, and is a no-op where the setting does not exist.
  Your global anndata settings are never changed.
• The CLI fixtures write through write_output too. The global setting stays
  off, so the CLI tests still exercise the CLI's own opt-in rather than a
  test-wide switch.

Coercing string columns to object was the other option on the issue. I did
not take it: it rewrites the dtype of user data in the output file, and the
opt-in matches what anndata 0.13 now does by default.

anndata < 0.11 cannot write a StringArray at all, with or without an opt-in.
Under pandas 2 those versions produce object columns, so #23 does not arise
there, and the new write tests skip on them with that reason.

#28 and #33: the provenance sha in a linked worktree

One root cause. A linked worktree's git dir holds HEAD; its branch refs and
packed-refs live in the shared repository dir named by the worktree's
commondir file. _resolve_sha looked only in the per-worktree dir and
returned None. It now looks there and then in the common dir, follows a
symbolic ref a bounded number of steps, and stays stdlib-only, so no git
binary is needed.

A None inside a checkout is now logged. _resolve warns once, at
resolution time, when it finds a git dir but cannot resolve the sha. That
answers #28's "failing to resolve a sha should be visible" without making
provenance able to fail a run.

On #28's second point: the packed-refs fallback was already correct in an
ordinary clone, and test_resolve_sha_from_packed_refs covered it. It was
wrong only in a worktree, where it looked in the per-worktree dir. Both cases
are now tested against real git in every shape (clone, branch worktree,
detached worktree, each loose and packed), asserting equality with
git rev-parse HEAD. Against base, exactly the branch-worktree rows and the
synthetic commondir cases fail, while clone and detached pass: the flip set
predicted before running.

A guard added after review: a wrong sha is worse than none. _find_git_dir
walks up to any enclosing .git. So a Kompot tree vendored or copied into an
unrelated repository, nested below that repository's root, whether in a plain
clone or a branch worktree, was stamped with that repository's HEAD. Before this PR the worktree variant returned
None; the worktree fix would have turned that into a confident wrong answer.
Provenance now accepts a git dir only when its .git entry sits at the package's
parent, which is where it sits in Kompot's own checkout. Residual: a
kompot/ package copied to the root of another repository still stamps that
repository's sha, because by position it cannot be told apart from Kompot's own
checkout. Two real-git tests plant the nested case, asserting that the walk-up alone would resolve the unrelated repo
and that provenance returns None. The reviewer's independent 14-scenario
real-git matrix gives the same result through the guard as without it.

The guard caught this PR's own test harness. The suite comparisons run each
commit from a git archive exported inside the clone, which makes it exactly
a Kompot tree inside another work tree. Under the guard,
test_source_checkout_resolves_sha_and_editable failed there, because its skip
precondition used the unguarded walk-up. In the earlier rounds it had passed
there by stamping the enclosing clone's sha: the defect itself, read as a pass.
The test now skips unless the package's parent holds a .git (eac37df). In CI's plain
checkout it runs and passes as before.

This also retires the one branch-only failure discussed in "Tests" below,
test_source_checkout_resolves_sha_and_editable in a linked worktree.

#30: the docs build now fails on markup that survives rendering

docs/source/markup_leak_check.py runs at Sphinx's build-finished and raises
if any rendered page contains an unrendered literal, role or directive. It
needs no rule per RST construct, and it runs on Read the Docs as well as
locally. tests/test_docs_markup_leak_check.py covers the three shapes quoted
on #30, the code-block-comment role, and a real sphinx-build that must fail
on nested markup and pass on clean markup. That last test was red with the hook
neutered.

This section's first version was the second false zero on #30, and the
counts below replace it. The earlier "zero" in "Verification" missed
notebooks/ because its scan was depth-1. The check as first pushed here then
missed a whole shape. Markdown [`x`](url) becomes nested RST that renders
as <code>x __: the inner literal is a tag, so stripping tags
before matching removed the evidence, and the check passed with 14 such
links live. One of them (LazyGeneCovariance, notebook 03) was added by this
PR. A fourth pattern now matches the leaked hyperlink tail, and a test plants
that exact rendering.

Measured with the final check on a build of the pre-fix base (ec4ce74):
25 leaked constructs in 5 pages. 9 are literals nested in bold, 2 are
links whose text contained a literal, and 14 are links whose whole text was a
literal. (The check reports 27, because those 2 links trip two patterns each.)
The first round fixed 11 of them and this revision fixes the other 14. The
build at the new head is clean with the hook active. By git grep against
19ee1a1, 5 of the 25 came from this PR's own notebook edits and 20 predate
it.

What the check covers now: an unrendered double-backtick literal, a role, a
directive, and an anonymous-hyperlink tail, on every page at any depth. What
it does not: markup whose leaked form carries no marker, such as a single *
of italic nested in bold, because a bare asterisk is too common in real text.

Boundary: it catches markup that should be markup. A single * of italic
nested in bold is not caught, because a bare asterisk is too common in real
text to flag.

#31: dask_num_workers is deprecated, not plumbed through

I checked what the knob could govern before choosing. The Dask covariance
tensor is returned lazily, and since #26 the Mahalanobis step materialises it
one gene at a time, so there is no pool of workers for any setting to
bound. Plumbing it into de() would have added a second knob that does
nothing. Instead, passing it now emits a FutureWarning that says so, and
the global dask.config.set and the "Configured Dask to use N workers" log
line are gone. A test asserts Dask's global config is byte-identical after a
disk-backed predict.

Review change: ~200 genes, not ~1 000, for the sample-variance pass

Per @katosh's review. Every place that suggests a gene budget for the second
pass now says ~200: the guide, cli.rst, simplified.rst, the README, the
kompot.de docstring example, the DE config template and notebook 03. Every
figure quoted beside a budget is priced for 200, not left at 1 000:

• per-gene memory: "Two hundred genes is roughly 75 GiB" at 5 000 landmarks;
• the timing table's column is 200 × the measured pinned per-gene time
  (3 s, 12 s, 49 s, 7 min);
• per-sample imputations: 0.09 GiB;
• the null-gene note: 200 + 2 000 is a 2 200-gene tensor;
• the n_landmarks table is re-run from dry_run=True at 200 genes:
  5.4 / 14.6 / 29.7 / 78.3 GiB in memory, and 2.5 / 2.6 / 2.9 / 3.8 GiB
  offloaded.

The re-run uses the same synthetic 20 000 × 20 000 input as the published plan
table, and as a positive control it reproduces that table's 200-gene row
(3.1 / 78.3 / 3.8) and 1 000-gene row (6.4 / 380.2 / 7.7) exactly. The
measurements elsewhere in this body that mention 1 000 genes describe a
specific stalled run and are left as measured. The measured case for 200 is
already in this body: at 200 genes with store_arrays_on_disk=True this branch
runs in 57.2 s at 1.80 GiB, where origin/main took 697.9 s.

Suite, set against set

Full suite re-run after the review fixes, at eac37df, the pushed head. It is
compared against base ec4ce74
(this branch before any of the new fixes, main merged in). Each arm is a clean
git archive of its commit, run on one Slurm node with BLAS pinned to 2
threads, in venvs pinned with uv --exclude-newer 2026-09-25. The comparison is
by test id, not by count:

arm	base	this PR	passed on base, not on PR	failing on PR
Python 3.11, anndata 0.10.9, pandas 2.3.3	2152 passed, 36 skipped	2179 passed, 40 skipped	1 (see below)	0
Python 3.12, anndata 0.12.6, pandas 3.0.6	49 failed, 2106 passed, 33 skipped	2183 passed, 36 skipped	1 (see below)	0
Python 3.12, anndata 0.13.4, pandas 3.0.6	2153 passed, 35 skipped	2183 passed, 36 skipped	1 (see below)	0

The one row per arm is test_source_checkout_resolves_sha_and_editable, which
now skips in these exported trees. That is the vendored-tree guard at work:
see "The guard caught this PR's own test harness" in #28 above. In a real
checkout it runs and passes.

31 test ids are new; none disappeared. On anndata 0.12.6, 47 of the 49 base
failures now pass. The other 2 now skip, and a failure turning into a skip
is the direction worth checking, so I checked it. Both tests call
pytest.skip("Requires palantir mocking") unconditionally, after their fixture
write. On base they died at that write before reaching the skip, and they skip
identically on base under anndata 0.13. The extra skips on anndata 0.10.9 are
the new write_output tests, which need the setting that version lacks.

CI runs a single anndata (0.13.4, whatever pip resolves), so the 0.12 arm
above is the only place #23 is exercised at all. Pinning an anndata 0.12 CI
leg would keep it that way; I have not added one, because changing the CI
matrix is the maintainer's call.

The story

Supplying sample_col to kompot.de() replaces the single shared posterior
covariance with one (n_landmarks, n_landmarks) matrix per gene, per
condition. Kompot summed those two tensors into a third dense array
before use, then added the shared covariance into it gene by gene with
__setitem__. So:

• store_arrays_on_disk=True did not keep anything out of memory — with
  np.memmap inputs the sum materialised the whole tensor in RAM, which is
  exactly the array the flag exists to avoid (store_arrays_on_disk=True does not reduce peak memory, and with dask nothing is written to disk #26);
• compute_mahalanobis_distances then ran cov = jnp.array(covariance) in its
  gene-specific branch and never read cov — a second full-size copy, for
  nothing (store_arrays_on_disk=True does not reduce peak memory, and with dask nothing is written to disk #26);
• the dry run, which the docs would tell you to run first, priced a run with
  zero null genes while the run it described would add 2 000 (dry run does not resolve null_genes="auto", under-estimating the plan by 2 000 genes #25);
• and none of the user-facing surfaces mentioned the cost at all, while
  cli.rst recommended the single most expensive configuration the package can
  express (docs: the cost of sample variance is not stated where users meet the decision #27).

What changed

kompot.utils.LazyGeneCovariance holds references to the per-condition
tensors and materialises one (n_landmarks, n_landmarks) matrix when the
Mahalanobis step asks for a gene. The dense sum and the __setitem__ loop are
gone, and so is the dead jnp.array(covariance).

A LazyGeneCovariance deliberately does not take the Dask delayed branch.
That branch submits every gene as a task at once, and a lazy view materialises
its own slices inside each task, so the two nest and peak memory scales with
concurrency rather than with one gene. This was not in the original diagnosis —
the first version of the fix left the Dask arm at 1 886 MiB where a run with no
sample variance sat at 972, and only re-measuring found it.

null_genes="auto" is resolved above the dry-run branch, and the estimator
now raises on a null_genes it cannot interpret rather than falling through to
zero — a str reaching that comparison is the failure signature.

The plan describes the code that runs: no third tensor, disk charged only on
the path that writes, and the per-gene working set charged explicitly.

Measured

Same script, same machine, same synthetic input, 19ee1a1 against this branch.
Peak anonymous memory from /proc/self/smaps_rollup, so page cache for a
memory-mapped file cannot flatter the reading. 1 500 cells, 4 donors.

150 genes, 600 landmarks (412 MiB per tensor):

run	0.8.0 anon	0.9.0 anon	0.9.0 Rss	written	wall
no sample variance	949 MiB	972 MiB	1 218 MiB	0 → 0	33.0 → 32.6 s
sample variance, in memory	2 401 MiB	1 973 MiB	2 238 MiB	0 → 0	45.7 → 46.0 s
store_arrays_on_disk, with dask	2 782 MiB	1 153 MiB	1 442 MiB	0 → 0	231.8 → 47.8 s
store_arrays_on_disk, no dask	2 284 MiB	1 139 MiB	2 224 MiB	824 → 824 MiB	48.2 → 47.7 s

Two instruments, because one is not enough. Anonymous counts private heap
pages; Rss counts those plus resident file-backed pages, which is exactly
where a memory map puts its data. As a share of the extra memory the in-memory
run costs: the dask path is 22% on both, and the dask-less path is 18%
on anon but 99% on Rss.

So the claim to rely on is the dask path, which wins on every instrument and
is the recommended one. The dask-less path converts unreclaimable anonymous
pages into reclaimable page cache — a real win for surviving a memory squeeze,
not a smaller resident footprint, and under a cgroup (Slurm --mem, a
container) page cache is charged to the same budget anyway. An earlier revision
of this PR reported only Anonymous and justified it as stopping page cache
"flattering the reading", which was backwards: choosing Anonymous is what
excludes those 823 MiB.

120 genes, 900 landmarks (742 MiB per tensor): 3 502 → 2 809 in memory,
4 341 → 1 325 with dask, 3 079 → 1 341 without.

Disk-backed storage is now cheaper than in-memory storage (1 153 against
1 973), where before it was dearer (2 782 against 2 401). In-memory cost
per gene drops from 3 × n_landmarks² × 8 to 2 × n_landmarks² × 8 bytes —
0.37 GiB rather than 0.56 GiB at the default 5 000 landmarks.

#25, measured: on a 4 000 × 500 input the default dry run reported 0.485 GiB
for a run the explicit null_genes=2000 plan priced at 2.156 GiB — 4.44×
optimistic. Now identical.

Results are unchanged. Distances differ only by floating-point summation
order: max relative difference 1.7e-12 across the in-memory, Dask and
memmap paths, which is the same order as the pre-existing difference between
those paths on 19ee1a1 alone (1.7e-13 absolute on values ~0.5).

One cost moved the other way

Without dask the tensors are memory-mapped, and a gene slice of a
C-contiguous (n_points, n_points, n_genes) map is strided, so reading one
gene at a time touches the whole file where the old code did one sequential
read into RAM. On a clean 900-landmark pair: 74 s against 66 s, about 12%
slower for 57% less memory, on the path that is not the recommended one.
Installing dask avoids it in both directions — that path is 4.9× faster
than 0.8.0 as well as lighter.

A first scale-2 run showed 646 s for that arm, but it was running under the
pytest suite at load 60+; the clean re-measurement is the 74 s above. Every
wall-clock figure here comes from a shared machine under load and should be
read as an order of magnitude. The memory figures are high-water marks of
allocation and are far less sensitive to that.

Read this first: what a skeptic pass changed after the initial review

w242sk returned check/second-pass with six findings. All six are accepted
and remediated in 03b4730; every figure below is a re-measurement of mine
with a positive control, not a copy of the skeptic's.

The largest finding invalidated a claim I had already rewritten this PR
around. An earlier revision said memory stops being the binding constraint
and compute takes over, on the strength of a timing table showing 8.5 s per
gene at n_landmarks=2000 and 47 hours for a transcriptome. That table
measured BLAS thread oversubscription on a loaded shared node, not Kompot.

It refuted itself, and I could have seen it without measuring anything: it
claimed 170.9 s for a 20-gene factorisation while a complete 20-gene
sample-variance run at the same n_landmarks finishes in ~48 s. A part
cannot exceed its whole.

Re-measured inside a real kompot.de run with OMP_NUM_THREADS=1:

n_landmarks	per gene (pinned)	per gene (unrestricted)
500	0.016 s	—
1 000	0.062 s	—
2 000	0.246 s	6.135 s
5 000	2.084 s	—

np.linalg.cholesky at n=2000 says the same thing from the other side:
21.3 GFLOPS pinned against 0.5 GFLOPS unrestricted. A throughput figure
that collapses as threads are added is contention, not work — that was the
tell, and it generalises well past this PR.

So the corrected picture, and the prescription rebuilt on it. 20 000 genes
is ~1.4 h at 2 000 landmarks and ~11.6 h at 5 000, not 47 h. Which means the
two-pass workflow rests on the memory argument, which is exact arithmetic
rather than a timing:

• held in memory, a whole transcriptome is 7 552 GiB — genuinely out of
  reach, and this is the hard limit;
• offloaded, it is 101 GiB and ~11.6 h — a long job, not a wall.

So ~200 genes is a budget rather than a barrier once you offload, justified
on value: pass 1 ranks the genes anyway, 200 costs about a hundredth on both
axes, and FDR is not calibrated for sample variance, so the other 19 800 buy no
testable calls. (This said ~1 000 until the maintainer's review; see
"Review change" above.) The docs say exactly that, and no longer imply a wall that is
not there.

There is a product finding inside this one, now a warning in the guide: the
per-gene loop is a single multi-threaded LAPACK call on a matrix that is small
relative to the core count, so on a busy shared node — which is the machine
most people run this on — it is ~25x slower than it needs to be. Pin
OMP_NUM_THREADS=1.

Tests

The first version of these tests had no potency, and the skeptic proved it

All 14 cases passed against a mutant that rebuilt the dense
(n_points, n_points, n_genes) sum in __init__ — #26's exact defect. They
asserted numerical equivalence and input validation, and an eagerly-materialising
implementation satisfies both, so they were never about allocation at all. Base
absence could not stand in either: the pre-existing suite passed on 19ee1a1,
which had that behaviour.

Two guards were added that assert the allocation property itself:

• a recording proxy that distinguishes [:, :, g] from any whole-array read
  (np.asarray, +, a slice spanning genes) — an eager build cannot avoid the
  second kind;
• a measured guard that walks a 62 MiB tensor and asserts anonymous memory
  grows by roughly one gene's matrix rather than by the tensor.

Both were verified red against the same mutation, the measured one reporting
"anonymous memory grew 67.4 MiB walking a 61.9 MiB tensor", while the other 14
stayed green on it. That last part is what makes this a demonstration rather
than an assertion: the new guards bite and the old ones still do not.

…and one of them was a false green in CI's configuration, which the skeptic caught

I asked the skeptic not to take the paragraph above on my word. It found that
test_peak_memory_stays_bounded_by_one_gene passed against the mutant whenever
the whole file ran, and failed only standalone — so the configuration I
verified in was not the configuration CI uses.

The cause was one line: the measurement window opened after the constructor.
An eager implementation materialises in __init__ — #26's actual shape — so
the 62 MiB allocation happened entirely outside the window, and what the guard
measured was incidental churn whose size depended on allocator state left by
sibling tests.

Reproduced before fixing (whole-file: only the proxy guard failed; -k peak_memory alone: the memory guard failed), then the window moved above the
constructor. Now red 3/3 whole-file against the same mutant and green 3/3
restored, with kompot/utils.py md5-verified against a pre-mutation backup.

The transferable rule: run a potency check in the configuration the suite
actually runs in, not the one that isolates the case. Isolation is what makes
such a check readable, and it is what made this one wrong.
kompot/utils.py was restored from a byte-checked backup afterwards, md5
verified.

The rest

tests/test_sample_variance_cost_regressions.py, 16 cases. Every expected
value is derived from the allocation (2 * L**2 * n_genes * 8,
3 * L**2 * 8), not copied from a run. Verified to fail on 19ee1a1: the five
that can run there fail, and the LazyGeneCovariance cases cannot even import.

One existing test changed, and it is the one predicted. The set of tests
expected to change was written down before the estimator was touched:
test_resource_estimation.py::test_dry_run_with_disk_storage, because it
asserted a bare len(disk_reqs) > 0 and disk is now charged only on the path
that writes. The failure-set diff after the change is exactly that test —
nothing else new, nothing gone. The prediction held in both directions.

It is re-derived, not loosened, and is strictly stronger than what it
replaced: the old assertion was an existence check that passed for the wrong
reason. The new one asserts an exact byte count from the allocation on the
no-dask path, zero on the dask path, and that the tensors are not charged
to memory in either. No assertion in this PR was widened — no == became a
>=, no exact figure became a range, no tolerance grew.

Suite result is reported as a set difference, not a count — a matching
total is not evidence in either direction.

Against the same tree before any of these changes, the failing set is
byte-identical: 26 failures in, 26 out, nothing added and nothing removed.
The one test the change did break, test_dry_run_with_disk_storage, is fixed
by the re-derivation above, so the branch ends with the failing set it started
with and 11 more passing tests.

Against a clean origin/main worktree with a package-for-package matched
venv (diffed with importlib.metadata, identical apart from pytest-cov):

--- ONLY ON BRANCH ---
tests/test_provenance.py::test_source_checkout_resolves_sha_and_editable
--- ONLY ON origin/main ---
(none)
--- COMMON: 25 ---

That one branch-only failure is not a regression, and it is measured rather
than argued. It is #28: _provenance.py cannot resolve a sha when HEAD is
a branch ref inside a linked worktree. The origin/main arm escaped it only
because git worktree add --detach leaves HEAD holding the sha literally. I
attached that same origin/main worktree to a branch and re-ran the single
test against origin/main's own code:

now on branch: tmp-provenance-check
HEAD file: ref: refs/heads/tmp-provenance-check
E   AssertionError: sha must resolve inside a work tree
FAILED tests/test_provenance.py::test_source_checkout_resolves_sha_and_editable

Detached again, it passes again. So the failure tracks the shape of the
checkout, not this branch. (The worktree was restored to detached 19ee1a1
and the temporary branch deleted.)

This is the only "only on branch" row, so it is the one a reviewer will stop
on, and it is worth saying what it turned into. It could have been explained
away with the #28 mechanism, which was already in hand — but an explanation
that flatters the branch is the one to distrust, so it was reproduced instead.
Doing that also confirms #28 from a direction #28 does not itself claim: the
issue reports that a sha fails to resolve in a linked worktree, and this shows
the resulting test failure appearing and disappearing purely as HEAD is
attached to and detached from a branch, on unmodified origin/main code.

Since fixed in this PR (#28 and #33, 04c8574). This row no longer
fails in a linked worktree on a branch.

The common failures are local-environment artefacts, not CI ones: 24 are
plotting tests that need scanpy, which neither local venv has, while
.github/workflows/tests.yml installs .[all] and therefore does have it. One
is the provenance test — see #28, and it fails here because these are linked
worktrees, which CI's plain checkout is not. One is a CLI diffusion-map test.
None are in code this branch touches, and I expect CI to be green where local
is not; if it is not, that is worth reading rather than assuming.

On reading that comparison: both arms ran back to back on the same shared
machine, with the 1-minute load average sampled throughout and ranging from
42 to 96 across the pair, which is high. Load
manufactures false reds, never false greens, so the greens are trustworthy
and any red outside the established set would be a statement about the machine
rather than about the branch. Both arms met comparable load by construction —
they ran sequentially, not concurrently — which matters because a spuriously
red test on the origin/main arm would read as "fixed by the branch", the
flattering direction and therefore the one to guard against.

Documentation

New docs/source/resource_planning.rst, the canonical explanation, linked from
index.rst, simplified.rst, anndata.rst, cli.rst, README.md and the
markdown of notebooks 02 and 03 rather than duplicated into them. It was
written from the post-fix measurement, not patched from a pre-fix draft.

The framing changed substantively as a result. Memory is no longer the reason
to restrict genes — disk offload answers that — so the two-pass workflow now
rests on the per-gene Cholesky, which is the one cost no storage mode touches.
Measured: 0.52 s per gene at 500 landmarks rising to 8.5 s at 2 000, against a
shared-covariance column that is flat in the gene count.

Also corrected:

• kompot.de's docstring, GPSettings.n_landmarks, GPSettings.batch_size,
  StorageSettings, SampleVarianceEstimator, kompot de --sample-col and
  the DE config template, so help() and --help carry the warning too.
• cli.rst's "Example: Complete Analysis" paired --sample-col Sample with
  --n-landmarks 5000 over every gene — copy-pasteable, and presented as the
  recommended complete form. Replaced with the restricted second pass through a
  config file, verified working end to end.
• StorageSettings.max_memory_ratio was documented as "Fraction of RAM
  before triggering disk storage". It triggers nothing: store_arrays_on_disk
  defaults to None, which resolves to disk_storage_dir is not None, and
  SampleVarianceEstimator.predict() hardcodes 0.8 for its own check, so a
  user's value never reached it. The config template's
  store_arrays_on_disk: null # (null = auto) implied the same non-existent
  automatic decision.

A finding worth more than its fix: the comparison snippet printed the inverse of the claim

The guide's with/without comparison — the snippet a reader is most likely to
copy, because it is the one that answers "what does sample variance cost me?" —
left null_genes at its "auto" default on both arms. After the #25 fix,
"auto" resolves to 2 000 without sample_col and to 0 with it. So the
two arms were silently priced over different gene counts, and the snippet
printed:

ratio = 0.5x

Sample variance looking half the price of no sample variance — the exact
inverse of the claim this entire PR exists to make, in the most copyable
artefact in the document. With null_genes=0 pinned on both arms it prints
5.8x.

Nothing in review would have caught this. The code was correct, the API was
used correctly, and the defect was entirely in a configuration that is
asymmetric by design and invisible at the call site. It was caught because
every published snippet is executed rather than read, which is the practice
that earned it and the reason to keep doing it.

It is also the second time on this branch that the "auto" sentinel produced a
confidently wrong number in the safe-looking direction; #25 is the first.

A claim withdrawn against my own measurement

I had written that the per-gene Cholesky is "roughly cubic in n_landmarks".
A Cholesky is O(n³) in flops, so it reads as obviously right, and it was in the
guide, in GPSettings.n_landmarks and in notebook 03's lever table before I
checked it against the table sitting directly above it.

500 → 2 000 landmarks is 4×, and 0.52 → 8.5 s per gene is 16×: an overall
exponent of 2.0, not 3. Worse, the per-interval ratios are 11.7×, 1.3×,
1.1× — exponents 3.55, 0.62, 0.30 — which is not a power law at all. It looks
like JAX warm-up at the first size followed by a flat region.

No exponent is claimed anywhere now. The four points are presented as four
measured points in the range they cover, with an explicit instruction not to
extrapolate and to time a single gene at your own landmark count if the number
matters.

The memory claim is untouched, and the reason is the part worth keeping.
2 × n_landmarks² × n_genes × 8 is derived from an allocation — it is the
shape of an array the source asks for, and it is exact. The cubic claim was
fitted to four timings. Those are different epistemic objects and should
never be stated in the same voice, however similar they look once they are both
sitting in a sentence as "quadratic in n_landmarks" and "cubic in
n_landmarks". One can be read off the code; the other needed evidence it did
not have. A plausible mechanism — a Cholesky is O(n³) in flops — is the
hardest kind of wrong claim to see, precisely because review confirms the
mechanism rather than the measurement.

A withdrawn claim that came back in its own remediation

6685622 withdrew "roughly cubic in n_landmarks" as unsupported. The
remediation commit then reintroduced it — "close to the cubic flop count a
Cholesky implies — 500 to 5 000 is 10x the landmarks for 130x the time" — a
sentence that refutes itself, since cubic is 1000x and not 130x.

The measured exponent over the four pinned points is 2.11, consistent
interval to interval (1.94, 2.00, 2.33). Cubic predicts 16 s per gene at 5 000
landmarks against 2.08 s measured, a 7.7x overshoot. It now says quadratic-ish
with the mechanism: the step carries O(n²) work that dominates at small n,
and single-threaded dpotrf throughput climbs 8.2 → 37.7 GFLOPS across the
range, cancelling much of the third power — with the note that the exponent
should drift toward 3 at larger n.

The irony is the content. The original withdrawal was right because the
unpinned data was not a power law at all, so no exponent was claimable. Pinning
the threads is what made a scaling statement possible for the first time — and
the statement made was the wrong one. A correction can restore the error it
corrected, once better evidence makes the claim sayable again.

Two further findings were consequences of the same fix rather than survivors of
it, and both are swept rather than patched: four sites (including de()'s
docstring, which ships in help()) claimed n_landmarks does not affect
per-gene compute, contradicting GPSettings.n_landmarks in the same
package; and the CHANGELOG — the file that becomes the release note — had
never been touched by the entire remediation, still carrying Anonymous as its
sole instrument and an unqualified "cheaper than in-memory".

A sweep that reported clean, and the instrument that made it lie

Three skeptic rounds; the third found that my post-fix sweep had missed three
sites still carrying retired claims — including, in
examples/03_sample_variance.ipynb cell 0, verbatim the sentence the README
fix replaced, in the notebook's opening cost warning.

The instrument is the finding. Every one of those sites is line-wrapped,
so a line-based grep cannot match it:

255:   matrix and its own Cholesky factorisation: about 0.37 GiB and several
256:   seconds per gene at ``--n-landmarks 5000``.

My search for 'several seconds per gene' returned zero against text that
plainly contains it, and the zero ran in the flattering direction — it
confirmed a remediation that was incomplete. The working form:

git show "${ref}:${path}" | tr '\n' ' ' | tr -s ' ' | grep -iF 'phrase'

Re-swept across all 219 tracked text files with twelve retired phrases,
including eight negative controls: zero hits. So the rest of the sweep was
complete and those three were the whole remainder.

The two files failed for different reasons, and the second is the one to
learn. docs/source/cli.rst was never opened by any remediation commit, so
the sweep never reached it. But examples/03_sample_variance.ipynb was opened
by the first remediation commit, to clear its stale outputs — and the
retired claim in its opening cell was left in place. A file being open in the
editor is not a file being swept.

Two related repairs came with it: cli.rst and kompot/cli/de.py had
disagreed about --sample-col (the --help text saying ~2 s per gene, its own
documentation page saying several seconds), now reconciled; and the CHANGELOG's
memory table had acquired an Rss column spliced from a different run than
its Anonymous column, so a reader could compute file-backed residency across
two runs and get 289 MiB against the 262 MiB the docs state. Rss is dropped
from the release note entirely and cross-referenced to the docs page — a
self-consistent single-instrument before/after is the right shape there, and it
cannot re-acquire the defect.

Verification

• sphinx-build -b html exits 0 with 24 warnings, an identical set to
  origin/main — no new ones. They are 22 pre-existing duplicate-object
  warnings from anndata.rst overlapping simplified.rst, one
  anndata.rst isn't included in any toctree, and one nbsphinx note about
  ipywidgets. (anndata.rst is orphaned on main too; out of scope here.)
• Rendered HTML scanned across every built page for raw reStructuredText
  leakage — unrendered literals, directives, unresolved roles: zero. Two
  real leaks were caught this way and fixed (nested inline markup in the de()
  docstring, a :ref: inside a code-block comment).
  Correction (2026-09-25): that zero was wrong, twice. The scan globbed
  {build}/*.html and never descended into notebooks/. Then the first
  version of the replacement check could not see code-in-link leaks. The
  measured total is 25 constructs, and all are fixed; see docs build does not catch nested inline RST markup (renders as literal backticks, sphinx exits 0) #30 above.
• Every code snippet published in the guide was executed against a
  synthetic AnnData — see the finding below for why that is not ceremony.
• All branches of the restructured assembly verified bitwise identical
  (maxdiff 0.000e+00) to the dense arithmetic they replace: both predictors,
  variance_predictor1 alone, variance_predictor2 alone, and a 2-D shared
  term, which the view refuses rather than broadcasting. Kept as tests.
• use_empirical_variance=True combined with sample_col verified separately,
  since the empirical diagonal is added to the per-gene matrix after the view
  assembles it: identical to the dense path (maxdiff 0.0), demonstrably
  changing the result rather than being silently dropped, and finite end to end.
• The two-pass field semantics in the guide were measured, not assumed: pass 2
  adds ..._mahalanobis_sample_var (NaN outside its gene subset), leaves pass
  1's ..._mahalanobis and ..._mean_lfc intact, rewrites the shared layers
  only for the genes it analysed, and adds two _std layers.
• Version is 0.9.0 in kompot/version.py, pyproject.toml and the changelog
  heading, and the installed import agrees.
• pyflakes over the whole package, diffed against origin/main as a set of
  (file, finding) pairs so line-number drift does not mask anything: identical,
  none added and none removed. It caught two dead imports this branch created
  (jax and tqdm in differential/differential_expression.py, orphaned when
  the eager assembly went), which are now gone. CI does not run flake8, so this
  was housekeeping rather than a gate.

Deliberately not in this PR

• CLI write path fails on anndata >= 0.11: allow_write_nullable_strings is False #23 (CLI write path on anndata >= 0.11) is older and unrelated to this
  story. Not folded in, so the omission is visible rather than accidental.
  Now bundled at the maintainer's request (2026-09-25); see CLI write path fails on anndata >= 0.11: allow_write_nullable_strings is False #23 above.
• kompot/_provenance.py cannot resolve a sha in a linked worktree checked
  out on a branch. It reads <gitdir>/refs/heads/<branch>, but a linked
  worktree's branch refs live in the main gitdir, so it returns None and
  stamps runs with a null sha. Pre-existing — confirmed by calling
  origin/main's own _resolve_sha against this worktree's gitdir — and it is
  the single head-only failure in the A/B, a property of the tree layout rather
  than of this branch. Filed as provenance: git sha silently resolves to None in a linked worktree (and packed-refs is a second silent None) #28, with the packed-refs variant that a worktree-only fix would leave live.
  Now fixed here, together with _provenance: sha does not resolve in a git worktree with a branch checked out #33; see "provenance: git sha silently resolves to None in a linked worktree (and packed-refs is a second silent None) #28 and _provenance: sha does not resolve in a git worktree with a branch checked out #33" above.
• The memmap slowdown above could be removed by storing the tensor gene-major,
  so a gene slice is contiguous. That changes the on-disk layout and is a
  bigger change than a minor version should carry; the recommended dask path
  does not have the problem.

An independent reproduction, and one number this PR does not yet carry

Added from a separate investigation that started from a stalled user run rather
than from this branch, and arrived at the same defect by a different route. The
run supplied sample_col over a restricted list of 1 000 genes at 1 000
landmarks with store_arrays_on_disk=True — the two-pass workflow this PR
documents, already followed — and did not finish, on a dataset whose pooled
first pass had scored every gene in 93 s.

What the graph actually does, as a closed form rather than an exponent

The __setitem__ loop this PR removes has an exact cost, and it is worth stating
as a count rather than as a time. Instrumenting the loop directly and reading
len(arr.__dask_graph__()) after each iteration, at 200 points and 120 genes:

 iteration       tasks
         1       2,891
        19      43,490
        59     276,390
       119     994,740

Least squares over all eleven points gives

tasks(g) = 61.5 g² + 1025.5 g + 1804      max relative residual 0.0%

An exact quadratic, not a fitted exponent — which matters here given the
withdrawal recorded above. A naive power-law fit over the same points returns
g^1.27 with a 61% residual, and the per-interval exponents climb 0.50 → 1.86,
because the linear term dominates until g is large. Reading the exponent off
successive intervals understates this one; the closed form does not. The
quadratic term is the __setitem__ chain, the linear term is the stack itself.

Through the real DifferentialExpression.compute_mahalanobis_distances,
intercepting the object handed downstream: 7 790 tasks at 20 genes against
314 360 at 80 — 40.4× for 4× the genes, where linear is 4× and quadratic 16×.
Different constants from the isolated series because the shape differs; the
point is that it is not a constant factor.

End to end through kompot.de() on 19ee1a1, synthetic data at the shapes
that run had (~9 500 cells, 6 samples per condition, 1 000 landmarks), one node,
one job, random_state pinned:

n_genes	store_arrays_on_disk=True	=False
25	49.9 s	47.9 s
50	53.4 s	50.2 s
100	95.0 s	53.4 s
200	697.9 s	63.8 s

No exponent is claimed from that table. The ratio per doubling grows rather than
holding, so it is not a power law, and these are four measured points in the
range they cover. About 45 s of each row is fixed cost.

What this branch does to those numbers

Same job, same node, same synthetic data, store_arrays_on_disk=True throughout —
the configuration that stalled:

n_genes	origin/main	this branch	peak, origin/main	peak, this branch
25	49.9 s	47.3 s	2.25 GiB	1.64 GiB
50	53.4 s	48.1 s	2.65 GiB	1.65 GiB
100	95.0 s	50.6 s	3.56 GiB	1.73 GiB
200	697.9 s	57.2 s	5.39 GiB	1.80 GiB

12.2x at 200 genes, and the memory column is the part worth reading: this
branch's peak moves 10% across an 8x increase in genes (1.64 → 1.80 GiB)
where origin/main's moves 140% (2.25 → 5.39 GiB). That is the lazy assembly
doing what it says — peak is one gene's matrix rather than the tensor — and it
is visible from the second row onward. Every checksum matches across both
columns, so this is the same arithmetic at a different cost.

On this branch, store_arrays_on_disk=True at 100 genes (50.6 s, 1.73 GiB) is
now faster and lighter than origin/main with the flag off (53.4 s, 3.84
GiB), which is the first configuration in which the flag is worth setting.

At 1 000 genes — the scale the stalled run was actually asking for — the trade
is visible rather than free, and it is worth stating in the same breath as the
12.2x above:

at 1 000 genes	wall	peak
origin/main, flag on	not run (projected hours)	—
origin/main, flag off	159.5 s	24.78 GiB
this branch, flag on	236.8 s	3.35 GiB
this branch, flag off	114.8 s	17.46 GiB

This branch improves both paths: flag off is 28% faster and 30% lighter than
origin/main's, and flag on goes from unrunnable to usable. Within the branch
the flag costs 2.06x the wall clock for 5.2x less memory at this scale — a
real trade rather than a free win, and the direction a disk mode should trade
in. Below ~200 genes the flag is faster and lighter, so the crossover sits
somewhere between; I did not locate it.

Results are unchanged across all of it. At 100 genes with random_state pinned,
the sample-variance Mahalanobis vector was saved from five arms — origin/main
flag on and off, this branch flag on and off, and a broadcast-add variant I
wrote and discarded — and compared element by element across all ten pairs:
max absolute difference 9.3e-13, max relative 2.8e-12. Not bit-identical,
and it should not be, since this branch sums its terms one gene at a time rather
than up front. That is a wider and looser check than the maxdiff 0.000e+00
reported above for the lazy assembly against the dense arithmetic it replaces —
separately built trees, end to end through kompot.de(), rather than one code
path on the same inputs — and the two agree about what they each measure.

This PR's fix, checked against that reproduction rather than against its prose

The guard from that work is added here as
tests/test_sample_variance_graph_scaling.py. It drives the real
compute_mahalanobis_distances and asserts on the object handed downstream, so
it exercises whatever the shipped code builds rather than a particular helper,
and it counts tasks rather than wall time so it cannot flake on a loaded CI
machine.

• On origin/main it fails: 40.4x for 4x the genes.
• On this branch it passes, unmodified.

It passes because LazyGeneCovariance carries no task graph at all, so the guard
finds nothing to measure and returns. The blowup is not bounded on this branch,
it is unrepresentable — which is a stronger property than the test was written to
check, and the test now says so explicitly rather than asserting a dask-specific
shape. Its first version asserted the representation and errored here with
AttributeError: 'LazyGeneCovariance' object has no attribute '__dask_graph__';
that was the probe's assumption, not this branch's behaviour.

Two things deliberately NOT added, both because measuring them said not to

• The closure capture in the retained Dask branch. compute_gene_mahalanobis
  is captured rather than passed, so Dask cannot see the covariance as a
  dependency. I had this filed as a second compounding mechanism and intended to
  fix it here. Measured against this branch with a bare 3-D Dask array — the
  shape tests/test_mahalanobis_approaches.py passes in — the graph is linear
  and so is the cost: 0.67 / 1.16 / 1.80 s at 25 / 50 / 100 genes, falling per
  gene. Against origin/main with the same linear graph, 0.24 / 0.88 / 2.18 s.
  Within a small constant factor of each other. The capture only ever mattered
  because it multiplied against a quadratic graph; with the graph linear it costs
  nothing measurable. Removing the __setitem__ chain was the whole fix, and I
  am withdrawing the "two compounding mechanisms" framing I published on store_arrays_on_disk=True does not reduce peak memory, and with dask nothing is written to disk #26.
• da.from_array(np.load(file_path), chunks="auto") at memory_utils.py:889.
  No mmap_mode, so the array is fully resident before Dask wraps it, and it
  looks like the answer to "why does the flag not reduce peak memory". It is not:
  DiskStorage.store_array and load_array are not reachable from the
  differential-expression entry points on either arm, on origin/main or on
  this branch: the Dask arm writes nothing and returns a lazy graph, the
  no-Dask arm writes a memmap and then sums two memmaps into a dense in-RAM
  array, and neither route calls them. (Every call site in this repo is a test,
  which supports that without establishing it — they are public module-level
  functions, so enumerating one tree cannot rule out an external caller.)
  Changing it would read as an answer to store_arrays_on_disk=True does not reduce peak memory, and with dask nothing is written to disk #26 while moving nothing the DE path
  reaches.

Unchanged by any of this

store_arrays_on_disk does not auto-enable under memory pressure — the
"critical" path at memory_utils.py:391 only escalates warning text, and the
flag resolves from disk_storage_dir is not None unless set explicitly. So for
anyone on a released version, where this branch has not landed, the mitigation is
simply not to pass store_arrays_on_disk=True: on the Dask path it creates a
storage directory, logs an expected size, and writes nothing to it.

Not merging — over to you.

[settylab-dotto-bot]

One measurement that is an argument about this PR's test suite rather than
about the fix, and that you have no way to see from here. It was produced by an
adversarial review of my own contribution, not by me.

Your 16 cost-regression tests survive the mutant that reintroduces the defect

Against an immutable export of 8bfcae04, with the expected flip set written
down before running:

The mutant. Revert the wrap site in
differential_expression.py — replace the LazyGeneCovariance(...) construction
with the pre-0.9.0 eager form, i.e. bind the summed tensor to the same name as
its input and add the shared covariance gene by gene:

gene_specific_covariance = combined_variance          # aliases its own input
for g in range(combined_variance.shape[2]):
    gene_specific_covariance[:, :, g] = (
        combined_variance[:, :, g] + combined_cov_to_add
    )

Result:

suite	against that mutant
tests/test_sample_variance_cost_regressions.py (16 tests)	16 passed
tests/test_sample_variance_graph_scaling.py (added in c4c7a3f)	fails, printing 40.4x for 4x the genes

Controls, so the mutant is not trivially detectable and the guard is not
trivially sensitive: an inert control edit survives the new guard, and an
eager Dask expression with a LINEAR graph also survives it — correctly, since
that variant does not reintroduce the graph blowup.

What that means for the merge

c4c7a3f is currently the only test in either tree that fails on an end-to-end
reintroduction of this defect. Your cost-regression suite is measuring real
things, but on this particular mutant it is not load-bearing, so merging the fix
without the guard would leave the regression uncovered by CI.

Scope, stated so it is not oversold: the guard asserts graph scaling, not
peak memory. The linear-graph eager variant passing is the demonstration — a
regression that restored the memory blowup while keeping the graph linear would
not be caught by it. Given this PR's framing, that boundary is worth knowing
before relying on it.

To reproduce: apply the mutant above to a checkout of this branch and run both
suites. The guard drives the real
DifferentialExpression.compute_mahalanobis_distances and asserts on the object
that method hands downstream, so it does not depend on which helper builds it.

[katosh]

@settylab-dotto-bot suggest only 200 instead of 1000 genes for the with sample var rerun.