Leiden 1.0.0 - #9
Open
ianfd wants to merge 20 commits into
Open
Conversation
added 20 commits
August 1, 2026 15:56
Spots and bins sit on a known grid, so neighbours are integer arithmetic — no spatial index, no distances, nothing approximate. Sort by (row, col), build a row directory, and each neighbour is a binary search inside one row. Offsets verified against 10x's Space Ranger docs and by replicating squidpy's algorithm: Visium is doubled coordinates, (row + col) constant parity, six neighbours at (row, col±2) and (row±1, col±1). Documented capture area (78 rows x 64 spots = 4992) is a test case. Two traps this defends against: - squidpy ignores the lattice and runs Euclidean kNN on pixel coordinates with a median*1.3 cutoff. Exact on contiguous tissue, but the cutoff is global, so at 50% lattice occupancy about half its edges are not lattice-adjacent. Working from the integers has no such failure mode. - raw doubled coordinates are anisotropic (in-row 2 units, diagonal sqrt 2), so metric methods silently drop the in-row neighbours and return a diagonal-only graph. visium_isometric_coords() rescales so all six sit at the pitch, and passing offset coordinates instead of doubled is now a hard error rather than a wrong graph. Builds CSR directly in two passes; an edge list first would cost ~700 MB at HD scale. 10.5M bins / 21.1M edges in 954ms, 591 MB. 18 tests: exact adjacency on small grids, closed-form edge counts, interior and hand-computed boundary degrees, parity, tissue holes, permutation invariance, symmetry, u32::MAX coordinates. Plus the load-bearing one — Leiden on a lattice produces zero disconnected communities, i.e. spatially contiguous domains.
Point-based assays (Xenium, MERFISH, CosMx, Slide-seq) position cells rather than slotting them into a grid, so neighbours have to be searched for. Both builders run over a uniform-grid index rather than a kd-tree. That started as a performance choice — tissue is close to uniformly dense, so bucketing beats traversal — but became a correctness one. kiddo 5.2.2 is unusable here on two counts, both found by tests in this commit: - KdTree panics outright past 32 points sharing a coordinate on one axis. Every Visium row is 64 spots at one y; every HD row is thousands. - ImmutableKdTree returns neighbour ids that do not match the coordinates it measured. Querying a collinear run, it reported distance 0 for a point five units away: the distances are right for the true neighbours, the ids are not. A graph built from it joins each cell to someone else's neighbours — wrong, and wrong in a way that still looks like a plausible clustering. The grid has no such failure modes. Duplicates and collinear runs are ordinary inputs, everything is exact, and it needs no dependency at all, so the spatial builders are unconditional rather than feature-gated. Radius is symmetric by construction. kNN is not, so Symmetry::Union (scanpy's convention, degree >= k) or Intersection (mutual only, degree <= k) reconciles it explicitly — emitting both directions and letting construction merge them would silently give mutual pairs twice the weight. Ties break on the lower index, so pixel-rounded coordinates still give a reproducible graph. Verified against O(n^2) brute force across seven degenerate point sets (collinear, all-identical, duplicated pairs, dense spike, far clumps) x three radii and three k. 500k cells: knn(6) 299ms, radius(30) 80ms.
Parameter-free neighbourhoods: no radius, no k, and the triangulation adapts to local density on its own. Verified exactly via Euler — a triangulation of n points with h on the convex hull has precisely 3n - 3 - h edges. A triangulation covers the convex hull, so it invents edges across ventricles, concave boundaries and between detached fragments. Joining two sides of a hole is exactly the error that yields a plausible-looking meaningless domain, so pruning is not optional. HullPruning::Adaptive drops an edge longer than `factor` times the local edge scale at *both* endpoints. Three details, each measured rather than assumed: - local scale is the 25th percentile of incident edge lengths, not the median. A node on a hole rim has several edges spanning the void, and those inflate its median until the artefact looks normal — the statistic hides what it is meant to expose. At the median, 6 of 42 fragment bridges survived at any factor. - compares against max of the two scales, not min. With min, a genuine boundary between a dense nest and loose stroma keeps only 19-40% of its edges; with max, 92-95%. - factor defaults to 2.0: the largest value that fully separates two detached fragments, the one unambiguous requirement. Costs ~5% of edges on uniform points, less on real tissue. It will not clear a hole 2-3 cells across; that needs ~1.5 and costs 16% everywhere, and at that width it is arguable whether the cells either side are neighbours anyway. Gabriel and relative-neighbourhood graphs are not the answer despite being parameter-free: their emptiness tests pass for a hole, so Gabriel leaves 19 of 42 fragment bridges and drops 21% of real edges, RNG drops 52%. Degenerate input (collinear, all-coincident) is a hard error rather than an edgeless graph, which would cluster into singletons and look like a result. 11 tests including Euler exactness, hole clearance, fragment separation, density-boundary preservation, permutation invariance, and contiguity of Leiden domains on the result.
The adaptive rule turns out to be Zahn's 1971 inconsistent-edge criterion: delete an edge "significantly larger than the average of nearby edge weights on both sides", which is the same both-endpoints shape used here. He reports "a factor of 2 usually means the separation is quite apparent", with worked examples from 1.3 to 2.6 — so the default of 2.0, arrived at by measurement, is also the canonical value. Records his limitation too, since it carries over unchanged: the criterion "doesn't detect one-way gradients however steep", so a smoothly ramping density gives it no discontinuity to find. Also notes what it would take to add Gabriel/RNG later: Lingas (1994) extracts either from an existing triangulation in O(n), and the exclusion region must be open, since the closed variant disconnects on ties and gridded coordinates are made of ties.
The existing tests prove each builder produces the graph it claims to. These ask the different question: does what comes out survive what real data does to you — cells missed by segmentation, centroids off by a fraction of a cell, a parameter picked slightly differently. Writing them turned up a limitation worth stating plainly. A bare spatial graph carries no expression, so the only thing it can recover is spatial structure. Physical separation it gets exactly: detached fragments and masked-apart tissue come back as connected components, with no resolution to choose, unchanged by 30% subsampling, by coordinate jitter, and across every parameter tried. Boundaries *inside* continuous tissue it does not get. There is nothing for the objective to anchor to, so it returns a tiling whose blob size the resolution sets. The tiling is spatially contiguous and looks exactly like a result, and it is arbitrary — resample and the boundaries move. My first attempt at these tests assumed otherwise and asserted domain recovery on planted density nests; it came back at ARI 0.09, which is the correct answer to the wrong question. So the file now asserts the limitation too: on homogeneous tissue the tiling must *fail* to reproduce under resampling. If that ever starts holding up, the assertion fires and says the docs need revisiting. 9 tests. The meaningful guarantee across resolutions is that a domain may subdivide an island but never spans two.
Written before implementing spatial domain detection so the target cannot drift to meet the result. Target: median ARI 0.46-0.52 across the 12 sections, against a non-spatial baseline of 0.38-0.43. The range is a real disagreement, not imprecision — BANKSY self-reports 0.518 from its own deposited artifacts, while the Genome Biology benchmark puts it near 0.46 running the same method on the same sections. Take the lower one. Records the independent medians for ten methods, the protocol rule that decides comparability (median ARI over resolutions yielding the correct cluster count, not best-of-sweep), and two traps: BANKSY's published ARIs are on smoothed labels via a SmoothLabels call visible only in its deposited code and never mentioned in the paper, and configuration can swamp method differences — DeepST scores 0.538 and 0.229 on the same section in two independent benchmarks. Also pins BANKSY's Visium parameters, including that lambda is 0.2 for Visium domain segmentation rather than the 0.8 used elsewhere, and that use_agf=TRUE means mean + gradient, so a mean-only implementation is not a faithful replication.
Fills in the self-reported vs independent comparison, which is the part that decides what to trust: GraphST 151673 self 0.635 indep 0.633/0.638 gap ~0.00 BayesSpace 151673 self 0.55 indep 0.550 gap 0.00 STAGATE 151676 self 0.60 indep 0.493 gap 0.11 BANKSY 12-slice self 0.518 indep 0.469 gap 0.05 Two of four survive independent replication. BayesSpace also disagrees between two independent benchmarks (0.550 vs ~0.40 on the same section), so independence alone does not settle it either. Verifies from code what no paper states: sample 2 (151669-151672) is annotated L3-L6 + WM, so 5 clusters, the other eight sections 7. Also that the benchmark drops unannotated spots before clustering rather than before scoring — GraphST's own tutorial does the opposite, letting them train the model and vote in spatial refinement, and hardcodes 7 clusters for every section including the 5-cluster ones. Expands the independent table to 15 methods and marks which are machine-readable from the benchmark's raw per-run files versus recovered from the figure. BANKSY is figure-derived and new in the published version, so 0.469 carries +/-0.01 and is not an exact published value.
…ence Augments each cell's expression with a summary of its neighbourhood so ordinary Leiden finds spatial domains. Verified differentially against banksy-py's own functions via committed fixtures, not against a reading of the paper — which matters, because four details are wrong if you transcribe from the paper, and three of them I got wrong first time. - The lambda budget is not split evenly across harmonics. Each successive harmonic gets half the weight of the one before, so at max_m=1 the mean takes 2/3 of lambda and the gradient 1/3, not half each. - Every block is z-scored per column before scaling, so lambda mixes standardised blocks. Without it the block with larger variance dominates regardless of lambda. - The gradient term subtracts its own neighbourhood mean before the phase sum. This lives in the reference's matrix builder, not its weight construction, so a matmul against the weights silently omits it — my first fixtures did exactly that and would have passed against a wrong reference. - The neighbourhood size differs per harmonic, and the two references disagree about how: R uses k_geom[m] (18 for both on DLPFC), Python multiplies internally as k*(m+1) (so 18 and 36). R produced the published numbers, so k is explicit per harmonic here rather than derived. Also documents, from source: the scaled_gaussian kernel has no factor of 2 in its exponent; R and Python disagree on the ranked kernel, which additionally crashes for m>0 in Python; and the reference stores azimuths as float32, which is why the fixtures agree to 1e-6 rather than tighter. 5 differential tests over 5 decay kernels, harmonics 0-2, lambda 0 to 1.
Completes the spatial layer. fuse() blends an expression graph with a spatial one over the same cells, alpha*a + (1-alpha)*b, as the union of both edge sets. Normalisation::TotalWeight scales each to unit total weight first, and is the default because without it the graph with heavier weights wins at every alpha and the parameter stops meaning anything — the same failure the BANKSY blocks avoid by z-scoring. There is a test that fusing a graph with a 1000x-scaled copy of itself gives back uniform weights when normalised, and is dominated when not. per_sample() runs a builder separately within each slice. Building globally and deleting cross-sample edges afterwards is not equivalent for kNN: a cell at a section's edge loses neighbours instead of taking k from its own section. There is a test asserting both halves of that — every cell keeps full degree under per_sample, and the filter-afterwards approach demonstrably starves some. It matters because sections are routinely stored overlapping in one coordinate frame, so cells from different slices can sit arbitrarily close. 8 tests including endpoint recovery, edge-set union, symmetry, sample labels being arbitrary, and single-cell samples.
Adds what the crate has actually been measured against, with the comparison named in each case rather than left implicit. Accuracy: 160 committed leidenalg fixtures (-0.09% mean modularity single-seed, +0.42% best-of-2), igraph (+0.03%), exhaustive brute force (93.5% exact optimum against igraph's 94.4-96.3%), banksy-py (elementwise to 1e-6), and PBMC3k where cluster counts track scanpy within one at every resolution and agreement with the authors' cell types is 0.8599 against scanpy's 0.8609. States plainly that the leidenalg mean is noise — both are stochastic heuristics and leidenalg's own two fixture seeds differ by up to 3.5% on hard instances — and that what the fixtures pin exactly is the quality function's definition, not the optimiser's luck. Performance is framed against scanpy/igraph as a C implementation called from Python, not as an interpreter-overhead comparison: 2-3x on PBMC3k, 1.33x geometric mean on synthetic graphs. Also records what is *not* measured: spatial domain detection has no comparative result yet, so nothing in the README should be read as a claim about it, and points at the pre-registered benchmark target.
The #[ignore]d tests were running nowhere. Three of the four are real coverage, not diagnostics: the Visium HD scale run is the only thing exercising the 10.5M-bin path, the Xenium run the only 500k-cell one, and the exhaustive optimality sweep over Bell(12) is the strongest evidence the optimiser actually reaches optima. They cost 11 seconds together, so there was no reason for the gap beyond nobody having added the job.
Leiden's inner loop is a scatter-gather over CSR, so at scale it should be
memory-latency bound and reordering nodes for locality should pay. Measured
rather than assumed, since the answer decides whether it is worth doing.
n morton bfs edge span (orig -> morton)
200k 1.07x 1.01x 5555 -> 314
1M 1.13x 1.09x 27778 -> 701
Locality improves ~40x and buys 13%. It does grow with size — at 200k the
membership array still fits in L3, so there is nothing to fix — but slowly, and
extrapolates to roughly 1.2x at 8M.
That ranks reordering below parallelism: refinement alone is worth ~1.4x for
comparable effort, and communities are independent during it so there is no
shared-state hazard. Modularity is identical across all three orderings, which
is the expected result of a pure relabelling and a decent check that the
experiment measured what it claimed to.
The README claimed 2-3x against scanpy on PBMC3k. That number included scanpy
converting its sparse matrix into an igraph object inside the timed call, so it
measured the pipeline, not the algorithm.
Adds examples/vs_igraph.rs and tools/igraph_bench.py, which write one graph both
sides read and time only the clustering call:
nodes ours igraph speedup
20k 0.030s 0.041s 1.38x
100k 0.166s 0.257s 1.55x
500k 1.605s 1.846s 1.15x
1M 3.457s 4.211s 1.22x
Geometric mean 1.32x, and modularity is identical to four decimals at every
size. The pre-existing 1.33x claim therefore survives a clean re-measurement;
the 2-3x one did not, and is now stated as what it is.
Worth noting the margin narrows above 100k. Both implementations end up
memory-latency bound on the same scatter-gather pattern, which is consistent
with the reordering experiment and says the remaining headroom is parallelism
rather than single-threaded tuning.
Runs N independent clusterings concurrently, each on its own graph. Nothing is
shared, so the only thing that can stop linear scaling is the memory system —
which makes it an upper bound on any parallel local-move, since that does share
state.
threads throughput vs 1
1 0.78/s 1.00x
2 1.44/s 1.85x
4 2.48/s 3.17x
8 4.01/s 5.14x
16 6.21/s 7.95x
Worth knowing because the igraph comparison showed the margin narrowing above
100k, which suggested a bandwidth wall that would make parallelism pointless.
It isn't one: latency-bound is not the same as bandwidth-saturated, and the
misses overlap fine across threads.
With 89% of a run parallelisable (level-0 local moving and refinement), that
ceiling puts the realistic target at 2.5-3.5x on 8 threads rather than the 9.3x
Amdahl alone suggests.
Refinement is constrained within each community of the working partition, so two
communities share no state at all. That was always true; the obstacle was that
the borrow checker could not see it, because everything went through one global
Partition.
Restructured so each community refines on local state:
- Objective gains delta_insert_view, taking the three aggregates a move actually
depends on (total weight, community strength, community weight). delta_insert
now reads those out of a Partition and delegates, so there is still exactly one
copy of each formula and the two paths cannot drift.
- refine_into() works a single community with local ids and local arrays, and
returns labels. The caller rebuilds one Partition from all the labels, so every
aggregate is exact rather than incrementally maintained across a split.
- rayon over communities when the feature is on.
One behavioural change: refinement randomness is now one draw per node, taken up
front in node order. Previously each select() pulled from a shared stream, so a
node's randomness depended on how many nodes before it happened to need one —
which couples communities that are otherwise independent and would make results
depend on how the work was divided. Clustering quality is unaffected; the
leidenalg differential fixtures still pass.
Measured on a quiet machine, clustering only:
nodes seq par igraph par/seq par/igraph
100k 0.150s 0.117s 0.199s 1.28x 1.70x
500k 1.113s 0.912s 1.469s 1.22x 1.61x
1M 2.568s 2.202s 3.410s 1.17x 1.55x
Note the earlier igraph figures in the README were measured while a background
job was running and were ~25% slow; the corrected single-threaded ratio is a
consistent 1.33x at every size, which is what the original geometric mean said.
Off by default because local moving is still sequential, so this is 1.2-1.3x
rather than the several-x a fully parallel implementation would give.
tests/thread_invariance.rs demands byte-identical labels across 1/2/4/8/16
threads over three objectives, two seeds and four graphs. That is the test the
pre-0.7 parallel path lacked: a staleness bug does not crash, it returns a
slightly different and entirely plausible clustering.
- neighborhood: 625 -> 252 lines. Dropped 21 println!s from library code, the
per-call rayon ThreadPool that only ever matched the default one, and ~180
lines duplicated between the kiddo and hnsw paths. The kd-tree path was also
writing *squared* distances into `distances` where hnsw wrote real ones;
both are real distances now. Connectivities are bit-identical either way.
- grouping: 255 -> 58 lines. `NetworkGrouping` abstracted exactly one type and
10 of its 12 methods were unreachable. Gone; `aggregate` takes `&VectorGrouping`.
- hnsw_rs `simdeez_f` is now x86-only. anndists imports simdeez::{avx2,sse2}
unguarded, so it could not build on aarch64 at all - `cargo test --all-features`
now works on Apple Silicon.
- shared test fixtures behind a `testdata` feature instead of two copies that
had already started to drift.
- dropped the unused num-traits dep, `collect_constrained`, a no-op
`free.clear()`, an `#[allow(unused_imports)]`, and an index loop that never
needed to be one.
- refine: one body for the parallel and sequential paths instead of two.
Long doc blocks nobody reads, gone. Kept the one line that carries the finding, dropped the essay around it. src comments 1226 -> ~900 with nothing factual lost. The Delaunay research - Zahn's factor, squidpy/Giotto defaults, why Gabriel and RNG don't help, the Visium loader gotchas - moved to docs/spatial-graphs.md, where it can be read once rather than sat in a header. Also: README says 1.0 now, matching Cargo.toml.
`parallel` gated refinement only, while the spatial builders threaded unconditionally - so a default build was parallel for graph construction and sequential for clustering, and there was no way to turn the former off. Now all eight rayon call sites go through `src/par.rs`, and `parallel` is default-on. Nothing regresses: gating spatial behind an off-by-default flag would have cost 2.6-3.2x on construction that users already get. 500k cells default single-threaded knn(6) 271ms 766ms radius(30) 72ms 233ms `default-features = false` now gives a build with no rayon in the tree at all, which wasn't possible before. Caveat: `knn` pulls in hnsw_rs, which threads on its own. Also dropped ndarray's `rayon` feature - nothing used its parallel iterators. The feature must change speed and nothing else, so `examples/feature_parity.rs` fingerprints clustering labels, modularity, and both spatial builders' full CSR including f32 weight bits. CI diffs a sequential run against a parallel one. Currently identical.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.