Leiden 1.0.0 - #9
Open
ianfd wants to merge 13 commits into
Open
Conversation
added 13 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.
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.