Changelog

All notable changes to this project are documented here. The format follows Keep a Changelog and the project adheres to Semantic Versioning.

Unreleased

Added

• Leiden passes: n_iterations on method="leiden" / "leiden-cpm", and clustering::leiden_iterated in Rust. Every pass after the first starts its first local move from the previous partition, as igraph's n_iterations does, and runs on its own seed derived from seed; the passes stop early once one returns its input, and a negative n_iterations repeats until that happens, capped at 100. The default is 1, which is the single pass of 1.0 and keeps its labels bit for bit: pinned digests of the core on the planted-partition graphs are unchanged. One pass stops at a local optimum of its last aggregate graph, and that is most of the detection gap reported below at s = 4, though not at s = 3. On the three-community SBM at s = 4 (medians of graph seeds 0/1/2) two passes move ARI from 0.657 to 0.897 at n = 3000 and from 0.563 to 0.870 at n = 30 000; until stable they reach 0.933 and 0.901, where python-igraph's n_iterations=-1 reads 0.910–0.923 on the same graphs. Modularity never falls from one pass to the next, and below the threshold the passes recover nothing. On the benchmark's real leaf graphs the gain is small: median ARI over one pass / until stable reads 0.728 / 0.731 on digits, 0.281 / 0.311 on MNIST and 0.046 / 0.046 on covtype. A new argument on clustering::leiden would have broken a module frozen under semver, and Method is frozen too, so a Rust caller who wants passes calls leiden_iterated; the numbers are in docs/MATH.md.
• CFTree::refit_leaves_beam — refit_leaves with the pass's own route widened, so the leaves are fitted to the partition a beam-wide route reproduces. refit_leaves is unchanged and is its beam = 1 case; a width argument on refit_leaves itself would have been a breaking change to a module frozen under semver, which is the promise working rather than an obstacle to route around.
• vmf is measured on real text, and the embedding row of the head table no longer leads with it. On 20-newsgroups (the bench's TF-IDF, reduced to 50 dimensions and L2-normalised, every row its own leaf, medians of three seeds) vmf reads ARI 0.073 against spherical-kmeans' 0.175, with 38 % of the rows in one cluster where the largest topic holds 5 %. An independent movMF EM collapses the same way and tying the concentration across clusters removes it, so the cause is the per-cluster κ — neither the implementation nor the summary. The measurement and when to prefer which head are in docs/USAGE.md, A concentration per cluster, and where it costs.
• Fuzz targets for the three entry points that take untrusted input — csr (validate_csr, then summarize_sparse on whatever it accepted), insert (CFTree::try_insert on arbitrary floats under every feature model, then refit_leaves_beam) and decode (a model's CBOR through CFTree::validate, then routed through and inserted into). They live in fuzz/, a crate of its own that is in neither the crates.io package nor the sdist, and run with cargo +nightly fuzz run -a <target>. What the first runs found is under Fixed.
• A closed-form target for kmeans, end to end through the tree. The population k-means cost per point of N(0, I_d) is d − 27/(8π) at k = 3 (three 120° sectors) and d − 2/π at k = 2 (two half-planes). Both partitions are Lloyd fixed points by symmetry; that the three sectors are the global optimum is supported numerically, not proven. A Rust test in tests/integration_api.rs and a numpy-only case in tests/test_python.py hold SSE/N within 0.01 of it at d = 2, 2·10⁵ rows, max_leaves 100 and 500, seeds 0/1/2: no labels and no reference implementation needed. The largest miss measured is 0.0035 (Rust) and 0.0037 (Python), against one standard error of 0.0024 / 0.0034 at k = 3 / 2. A head that seeds k − 1 centres (1.361 at k = 3) and one that never moves its centres (0.976) both fail it.
• A planted-partition check of the Leiden core. On the symmetric sparse three-community stochastic block model, weak recovery is possible exactly when s = (a − b)²/(3(a + 2b)) > 1 (OpenAI, The exact reconstruction threshold for the three-state symmetric channel, 2026 preprint, Corollary 1.2). With the communities in contiguous node-id blocks, mean degree 6 and n = 3000, the in-module test of community.rs holds two controls over seeds 0/1/2: at s = 0.6 an ARI above 0.005 could only come from the node ids (measured 0.0015–0.0021; at most 0.0037 over 40 graphs), and at s = 5 the blocks must come back above 0.75 (measured 0.847–0.886). A local move seeded from id blocks of 100 nodes fails the first. Measured along s at n = 30 000 (medians of three graphs), the modularity core reads ARI 0.017 at s = 3, 0.563 at s = 4 and 0.798 at s = 5: it detects from about four times the threshold signal, where the information-theoretic limit is s = 1. The test reaches the core directly; the leiden head builds its own kNN graph, so Python cannot feed it an SBM.
• What the heads can and cannot guarantee, a section of docs/MATH.md that names each objective exactly. The 2-way RatioCut is N times uniform sparsest cut, which is NP-hard to approximate within any constant, so the cut unnormalised spectral clustering relaxes has no polynomial constant-factor algorithm unless P = NP (the spectral head relaxes NCut, which the result does not cover). Metric k-median has the exact threshold 1 + 2/e, which neither dc-median (exact in its ultrametric) nor kmedoids (squared Euclidean, not a metric) is subject to. The kernel clustering hardness below (8π/9)(1 − 1/k) is about an objective no head optimises. Maximum modularity cannot be approximated within any positive factor unless P = NP, and the Leiden core detects from about four times the three-community threshold. The OpenAI sources are cited by title, date and URL as 2026 preprints, and each cited result says whether its source machine-checks it in Lean: the propeller bound, the k-median threshold and the Unique Games theorem are checked; the sparsest-cut and kernel-clustering hardness and the three-state impossibility direction are not.

Changed

• leaf_refit now routes at route_beam instead of always greedily. A Lloyd step is only a Lloyd step against the rule it is read back with. Setting both knobs used to hand back leaves fitted to a greedy partition while every query the estimator performs — predict, predict_proba, assign_microclusters, outlier_scores — used the wide one; that is not a cheaper answer, it is a mismatched one. Two consequences for callers who set both: the tree itself moves, so every head sees it including kmeans and the mixtures, which are otherwise width-blind; and the pass costs a linear factor in the width on its route, on top of the 1.8–2.5× a plain fit one narrow pass already costs. At the default route_beam=1 nothing changes — beam <= 1 is the plain descent, and the leaf statistics are bit-identical to 1.0. The free fit_predict / fit_predict_sparse engine functions take no width and keep refitting greedily. Measured against a build with the change reverted, over six cells (digits\@90 and mnist-10k\@200 × kmeans/ward/gmm, ARI the median of seeds 0/1/2): three cells up, three down, mean +0.013, and the split follows the argument rather than cutting across it — the two ward cells, the only head that assigns by microcluster, gain +0.061 and +0.052, while the two heads that label from their own centres move −0.030/+0.034 and −0.003/−0.038. Worst case first: on mnist-10k ward the greedy pass under a width of 8 was destroying 0.065 ARI, more than the width had bought. The table is in bench/RESULTS.md.

Fixed

• memory_budget_mb no longer under-sizes the divisor it uses, so the budget is met instead of overrun. The per-leaf cost is now measured rather than derived from the struct definitions: the slope of process RSS against the realised leaf count, at four leaf counts per cell, for each feature × dim in {2, 20, 54, 784} (r² 0.978–1.0000). The old formula under-predicted every cell, from 1.15× to 19.5× — so a 512 MB budget bought a tree of 590 MB at best and 10 GB at worst, and it never erred the safe way. Two causes. A resident leaf is more than its declared arrays: its cluster feature is also folded into every ancestor node's feature, a rebuild allocates a second tree whose small chunks glibc keeps in its arena rather than returning to the kernel, and every Vec costs a header plus malloc rounding. And feature="fd" was costed as if it were "diagonal", ignoring the Frequent-Directions sketch entirely — an FdSketch holds min(32, dim) sketch rows of dim besides the mean, which is the whole of the 19.5× at dim = 784. The new constants sit 1.02–1.37× above every measured cell, deliberately: a divisor that under-predicts blows the budget silently, one that over-predicts hands back a slightly smaller tree than the caller could have afforded. A given memory_budget_mb therefore resolves to fewer leaves than in 1.0, which changes labels for callers who set it — max_leaves is unaffected.
• linalg::orthonormalize_rows (Rust) documented the opposite of what it does with a dependent row. It said such a row is left at zero. Its residual is rounding noise rather than zero, so it comes back a unit vector orthogonal to the rows above it — which is what every caller in the crate needs, since each is a subspace iteration in which a zero row is a fixed point. The behaviour is unchanged and now asserted; only a residual below 1e-150 is zeroed.
• A rebuild could silently drop rows. Six finite rows — ±1e200, ±3e200, 0 and 5 — fitted at max_leaves=2 kept three. When every merge candidate in a leaf cost +inf (a centroid gap that squares past f64::MAX) or NaN, the entry was paired with itself, merged into itself and deleted. A candidate is now always two entries and a NaN cost ranks last. A gap that overflowed no longer grows the threshold either: it carries no scale, growing to it would absorb every later row, and +inf is the one threshold load refuses.
• fit_predict_sparse panicked on a finite value past ~1.3e154 once the leader budget was spent, with index out of bounds: … the index is 18446744073709551615. The value squares to infinity, every later row sits at +inf from its leader, and the nearest-leader scan — < against an infinite start — found none, so the forced branch indexed usize::MAX. The scan now returns an existing leader whenever one exists, and the capped rescan finds an open leader at any distance, as its documentation always said.
• Betula(feature="fd") panicked on one float32 value past ~1.8e19, under every head: the square overflowed the sketch's Gram matrix, and the NaN eigenvalues reached a partial_cmp().unwrap() sort. The shrink already read every eigenvalue as max(λ, 0); the median it shrinks by now does too.
• load accepted four tree shapes that the next call hung or panicked on — a child that does not name its parent (an insert climbed a parent loop forever), a root that names one, a reachable node with no children, and an entry listed twice (the next rebuild merged it into itself and emptied a leaf). Each is now refused with the index at fault. Every tree this library builds satisfies all four, so no file it wrote stops loading.

The first three change nothing for data whose squares stay finite. Each fix has a unit test that fails without it, and every one came out of the fuzz targets or out of reading the code a finding pointed at. - leaf_cap=1 never held a tree to max_leaves. A rebuild merges two entries that share a leaf node, and a leaf holding one entry has no partner, so every rebuild merged nothing and every insert past the budget rebuilt again: Betula(leaf_cap=1, max_leaves=50) kept all 2400 rows of four blobs as leaves. When no leaf node holds two entries, the rebuild now pairs across the leaves of one internal node and then rebuilds the node structure, since those merges empty leaves. At leaf_cap >= 2 a rebuild always finds a pair inside a leaf, so no such fit changes. - branching below 2 and leaf_cap=0 are refused by Betula, fit_predict, BregmanBetula and WindowStream, with a ValueError naming both values. They used to fit — leaf_cap=0 without honouring max_leaves — and a Betula fitted with them saved a model that load then refused, as it always has. - A finite row whose squares a fit would overflow is refused by every method that takes rows, dense or sparse, with a ValueError naming the row. The bound is dtype.max / 2⁵⁶ on the squared norm — about 4.7e21 for float32, 2.5e291 for float64 — which keeps every sum of squared distances a fit takes inside the type for any row count up to 2⁵³. Past it a row went through the finiteness check and then: on a float64 row at 1e200, kmeans indexed out of bounds (a PanicException, through Betula.fit on sparse input as well) and normalize=True scaled the row to zero; 3000 float32 rows spread over ±4.5e18 left the radii of every non-directional head non-finite and made kmedoids panic. A float64 chunk that a float32 model casts down is checked after the cast, where 1e39 is inf. Needing no row count has a price: a small float32 fit past the bound that happened to stay in range is refused too — rescale, or pass float64. A bound on rows × largest square was the alternative, and it would have made the domain of partial_fit depend on the rows that came before it. The check is at the Python boundary; CFTree::try_insert in Rust takes any finite row, as before. How large a row can be. - The HDBSCAN blobs-contest numbers are re-measured on the tree 1.0.0 shipped, and the automatic min_samples reads 0.822 at 8 000 leaves, not the published 0.896. The rule was measured on the tree before 1.0.0's cost-ranked rebuild, and the rebuild's own re-measurement covered the quality tables but not the contest section of bench/RESULTS.md, docs/USAGE.md or the AUTO_MIN_SAMPLES_LEAVES doc, so all three quoted a tree that was never released. The rebuild (7e5471f) is the whole cause: built as a wheel beside its parent, the parent reproduces the grid the numbers came from to the fourth decimal and the commit reproduces today's build on every cell. On 1.0.0 the rule reads 0.821 at 2 000 leaves and 0.822 at 8 000 (medians of seeds 0/1/2; the fixed ten 0.612 and 0.717), against fast_hdbscan's 0.910. The drop is all seed 0, which loses one cluster at each budget (0.790 → 0.561 at 2 000 leaves, 0.896 → 0.799 at 8 000; the 2 000-leaf median is another seed's and did not move). On the six quality fixtures its gain on blobs is gone — the fixed ten no longer fails there, 0.423 against 0.429 — and on varied it grew, 0.568 → 0.839. Over the wider grid it still wins or ties in 16 of 18 cells, but the two losses are no longer between two failures: on two overlapping blobs at N = 100 000 the fixed ten reads 0.484 at 200 leaves and 0.452 at 2 000 against the rule's 0.000 and 0.337. Below the leaf mass the head is single linkage on either tree, and it now fragments where it chained: the contest's min_samples=10 call finds 24 clusters at 100 000 points and 88 at 500 000 instead of 3, at ARI 0.480 and 0.490 against the published 0.478. The constant is unchanged — 5 to 40 leaves are still inside the seed spread at both budgets — and no code changed. The contest table's time and RSS columns are the 2026-08-25 run, and now say so. bench/RESULTS.md.

1.0.0 — 2026-09-10

1.0 is a promise about compatibility, not a claim about features. Three things are frozen from this release on, and the rest of the crate is explicitly not:

• The Python package is the supported surface. Every keyword of Betula, every free function and every fitted attribute is covered by semantic versioning. python/betula_cluster/__init__.pyi pins the signatures and mypy.stubtest checks the compiled module against it on every run of the gate, so a wrapper signature cannot drift without failing.
• The Rust crate's contract is thirteen modules — tree, feature, distance, bregman, model, clustering, types, sparse, order, coreset, stream, window, validity. Everything else is #[doc(hidden)]: still reachable, because this repository's own benchmarks call into it, and free to change shape in any release. cargo semver-checks now runs against the last version published on crates.io on every pull request.
• A model file written by 1.x loads in every later 1.y. save writes a gzip-framed, version-tagged CBOR document; load accepts one schema version, refuses every other by name, and validates the tree it decodes instead of trusting it.

What the promise costs is stated rather than hidden. cargo semver-checks reports six breaking classes against 0.8.0: three from curating the surface (enum_now_doc_hidden, function_now_doc_hidden, struct_now_doc_hidden) and three from signatures that gained a parameter or moved to a try_ form (function_parameter_count_changed, method_parameter_count_changed, inherent_method_missing). On top of those, an f32 tree past 2²⁴ rows now computes a different — correct — weight and variance, where it used to stop counting. All of it is below, and this is the release to take it in.

Explicitly not part of 1.0: being fastest, having every head, or matching a specialist library on its own contest. README.md names the three cases where a different tool is the right answer.

Added

• A "choosing a head" table in the README, and the three cases where the answer is a different library. docs/USAGE.md has documented every head for several releases and its twenty-row table answers "which one" — but only once you are already reading it. The README said which datasets betula-cluster is for and never which of the twenty-seven method values to type, so the first question a new reader has was answered two clicks away. It now carries a five-family table (centroid / probabilistic / density / graph / directional) naming a head to start from and what to reach for next, above a link to the full one.

The exclusions are written in the same place rather than left implicit, each with the measurement behind it. N in the low thousands: the tree stops compressing — on digits (1 797 points) it holds 1 797 leaves and rebuilds zero times, at max_leaves of both 2 000 and 4 000, one leaf per point — which is why that row is a tie (0.467 against scikit-learn's 0.468) rather than a win. Density structure that overlaps: fast-hdbscan recovers the 100 k blob fixture at ARI 0.910 against our 0.478, and being 9× faster on half the memory is the trade, not a rebuttal. And a mixture whose likelihood has to be a raw-point likelihood: responsibilities are tied within a leaf and every ln p is computed from (n, μ, S), so a per-point density, a BIC against a model fitted to raw points, or a likelihood-ratio test are not what comes out.

• n_init — the k-means restart count is a parameter. Lloyd converges to a local optimum of the inertia, so the k-means heads have always kept the best of four k-means++ draws; the count was a constant. On MNIST-20k it is worth choosing: 25 restarts read ARI 0.3303 against the four-restart default's 0.3069 and sklearn.cluster.KMeans's 0.3244, in half scikit-learn's wall clock and 3.0× the default fit's (median of seeds 0/1/2, the published harness). It is not a free win and the docs say so — on digits the same 25 restarts score below 10, and on covtype-20k the column has no direction, because the draw is selected by inertia and inertia is not ARI. The default is unchanged, so no existing fit moves.

Only the heads whose labels come from an inertia-selected k-means take it (kmeans, spherical-kmeans, spectral, and the COP-KMeans behind fit_constrained). The EM heads select their restarts by likelihood, and the same sweep over the GMM restart count was measured non-monotone on all three datasets, so n_init on a head that would ignore it raises a ValueError instead of passing silently. Rust: Model::fit and the kmeans_auto / spectral entry points take the count (0 = the default), and KMEANS_N_INIT is public.

• The ELKI tree gap has a mechanism: how the final leaf set is chosen. Q5 left our default D0/D0 geometry 7.6–8.5 % behind BetulaLeafPreClustering in WCSS at a matched leaf count, with three candidates on file (split rule, rebuild threshold, absorption tie-break). It is none of those: it is how CFTree::rebuild chooses the final leaf set — it merges the closest pairs, and only pairs inside one leaf node. Priced exactly in cluster-feature arithmetic — a merge-only compaction is a closed form, so an offline globally-cheapest-pair merge is a bound — a summary of exactly the same size built from a 2× finer tree reads 0.72× our WCSS on digits and 0.68× on covtype-50k, i.e. 0.78× and 0.73× of ELKI's. The 8 % is the small end of what the lever is worth. Not a free fix: the finer build holds 2× the leaves and the global merge is O(m²), which is the trade the current policy was chosen under. Written up in research/RESULTS-estep.md with the numbers and the shape of a fix that keeps both bounds. E13 then separated the two halves of that policy: the ranking is 21 of the 28 points and the locality is 7, so the fix in Changed above needs neither the finer build nor the global merge.

• absorb="chi2" is documented as an operating point, with the measurement that says it is not a general win. T22 turned up the gate reading 0.1199 with 46 leaves on covtype where the six radius criteria read 0.086 with ~3 600, which looked like a lever nobody had written down. Swept over three datasets × three heads × three budgets (medians of seeds 0/1/2) it is not: the gate wins on covtype-ward (+0.034) and on tight budgets (MNIST-10k ward at 500 leaves, 0.1601 → 0.2834), loses on digits-ward (0.6428 → 0.4456) and on MNIST at a generous budget (0.4275 → 0.3209), and is mixed on the k-means and GMM heads. What generalises is that it chooses its own resolution — 46 leaves on covtype and 920 on MNIST whether you ask for 500 or 4 000 — which makes it a cheap operating point rather than a better criterion: the MNIST cell above is 17× the wall clock (7.38 s → 0.42 s) for −0.107 ARI. docs/USAGE.md now carries the table.

• The covtype / ward benchmark loss has a mechanism on record, and it is not the one that was filed. The cell (betula-ward 0.086 against sklearn-birch 0.131) was attributed to mass imbalance; Q1 measured that absent, and this separates the rest. Not the leaf budget — 500 to 16 000 leaves reads 0.0857–0.0900, flat over a 32× range. Not the leaf feature — spherical and diagonal agree to four decimals, as the ward linkage's inputs say they must. Not the absorption criterion within the radius family — six of them read 0.0859–0.0864 — though the mass-invariant chi2 gate reads 0.1199 with 46 leaves, the rival's score at a 78× smaller summary. It is the fixture: bench/comprehensive.py standardises the 20 000-row subsample, bench/_worker.py standardises all 581 012 rows and then subsamples, and the two draws differ as well. Crossed, ward reads 0.0861 / 0.1176 / 0.1059 / 0.1416 while sklearn-birch sits at 0.1265–0.1306 throughout — the sensitivity is ours, and the 44 binary one-hot columns carry it. The scoreboard cell stands (the rival is measured on the same rows), with the new mechanism written next to it in bench/RESULTS.md.

• The spectral head's restart count is measured, and the default stays 4. n_init became reachable on method="spectral" in 0.8.0's successor without ever being swept there. It is the one head where a restart is nearly free — its k-means runs on a k-dimensional embedding of at most max_leaves rows — and the measurement says so: a hundred restarts cost 0.20 s → 0.28 s on digits and 0.39 s → 0.44 s on covtype-20k, putting the default four at 0.5–2 % of the fit. It also says the restarts buy nothing, which is why the default does not move: moons 0.9999 and circles 1.0000 are identical to four decimals at 4, 10, 25 and 100, digits moves +0.0009, and covtype-20k's −0.0103 sits well inside its own 0.0657–0.0992 seed spread. Documented in docs/USAGE.md next to the k-means table, where the curve is not flat.

• betula --n-init N (CLI). The binary passed the 0 sentinel, so --method kmeans was stuck on four restarts while both Python entry points could choose. It carries the same head rule as the library: --n-init on gmm, gmm-full or ward is an error rather than a silently ignored flag, and --n-init 0 is refused because zero is the engine's "take the default" sentinel and not a fit anyone means to ask for.

• simplified_silhouette — the silhouette a summary can actually carry. validity() now reports a fourth index: Hruschka's simplified silhouette, mass-weighted over leaves, measuring to the cluster centroid rather than to the members. The classical silhouette needs Σ‖x − y‖, which is degree 1 in the norm and provably not a function of a cluster feature (a Z3 witness: {−1, −7/12, 0, 1} and {−4/3, 0, 0, 3/4} share n, Σx, Σx² and differ in Σ|xᵢ − xⱼ| by more than 0.1), so this is a declared surrogate, not an approximation that tightens.

The two losses are measured separately and are not the same size. The summary costs 0.0004 on covtype-20k at 5.6 points per leaf, at most 0.013 over the k-grid on 784-dimensional MNIST, and exactly nothing on a tree with one leaf per point, where a unit test pins it to the point-level definition. The surrogate class shifts the level a long way — 0.4393 against scikit-learn's 0.1420 on digits — so the number ranks configurations of one dataset and must not be compared to a published silhouette_score. Rank agreement with a 10 000-point sampled classical silhouette is τ = +0.71 on covtype, +0.62 on digits and −0.24 on MNIST, where the point-level surrogate disagrees harder still (−0.43): the disagreement is the surrogate, not the summary.

tune() can now select on it — objective="simplified_silhouette" or "medoid_silhouette" — scoring off the leaves at 39 ms per trial on MNIST-20k against 1456 ms for the sampled point silhouette, with no sampling seed in the answer. Rust: validity::simplified_silhouette is public.

• balance="auto" — the mass cap as a decision the data makes. balance bounds how much of the total mass one leaf may hold, and it has a known domain: over 27 cells (digits / covtype-20k / mnist-10k × kmeans/ward/gmm × three budgets) a fixed cap gains +0.08 to +0.25 ARI on all six cells where the heaviest leaf holds more than half the mass, and stays inside seed noise on the twelve below 0.1. That statistic is now readable by the estimator itself: "auto" summarises once, reads the share off the finished tree and, only if it passes 0.5, summarises the same rows again with the cap on from the first point. Re-measured over the same 27 cells, it is bit-identical to balance=None on all 21 non-firing cells and lands on the balance=4.0 answer on all six that fire (+0.0785 to +0.2514). The price is one extra pass, paid only where it fires.

partial_fit cannot summarise the same rows twice, so a stream gets a weaker mechanism instead: the tree watches its own mass distribution as it builds and arms the cap when a leaf passes half the mass. A cluster feature does not split back into points, so a leaf that has already absorbed the core keeps it — measured on mnist-10k at 250 leaves the streaming arming recovers +0.006 to +0.037 where the two-pass recovers +0.167 to +0.251. Both are documented in docs/USAGE.md; the default is still None.

Not everything the parameter was hoped to fix is a mass problem: the published covtype/ward cell (ARI 0.086 against sklearn-birch's 0.131) has a heaviest-leaf share of 0.041, and a fixed cap moves it by −0.0001. "auto" correctly does nothing there.

• CFTree::top1_mass (Rust) — the heaviest leaf's share of the mass, the diagnostic the docs already told readers to compute by hand.

• A float32 tree now warns before its weights stop counting. An f32 leaf saturates at 2²⁴ = 16 777 216 points: binary32 spaces its values 2 apart from there upward, so a unit-weight row leaves the weight exactly where it was. Measured, 16 781 312 identical float32 rows into one leaf report a weight of 16 777 216 against an exact float64 control. The mean and the scatter do not stop with it — the mean becomes an exponential moving average of span 2²⁴ and the variance inflates in proportion to the rows the weight dropped (1.0105 against 1.0000 after 200 000 rows past the ceiling) — so the summary degrades silently rather than erroring. The estimator now raises a UserWarning naming the heaviest leaf's mass when it passes 2²³, one doubling short of the ceiling and while every summary is still exact, and docs/USAGE.md documents the limit and the three ways out. This is a warning, not a fix: the accumulators are still f32, and widening them is a persisted-format change.

• assign::AssignPlan — pruned exact nearest-centroid assignment (Rust). A squared Euclidean distance is a sum of non-negative terms, so any prefix of it is a lower bound on the whole, and a candidate whose prefix already exceeds the best distance so far cannot win. AssignPlan walks the dimensions in stages and drops candidates as soon as that happens. It returns the argmin a full scan returns, ties broken by the lower index; there is no error budget and no parameter to tune.

Two levers, both measured. A hint (the row's label from a previous iteration) starts the threshold at that candidate's distance instead of infinity, which takes the dimensions read from 45–60 % down to about 20 %. A dimension order by descending between-centroid variance prunes soonest — on mnist at k=10, 42.7 % of the dimensions read against 66.6 % in the natural order.

Fewer reads is not automatically less time and the module documentation says where it is not. cargo bench --bench assign on a Ryzen 7 5800HS, single-threaded, hinted, against a full scan: 0.3–0.8× below d = 128 (the branchless AVX2 row kernel wins outright), 1.8–2.1× at d = 784–1024 for k ≥ 100, 8.2× at d = 1024, k = 1000, and 3.0–3.9× when the clusters are well separated. The dimension order is a k decision rather than a d one: below about k = 100 the per-point gather costs more than the extra pruning buys, and identity_order is faster. Nothing in the library calls it yet — this is the kernel, and each call site gets its own measurement before it is wired in.

The statistical alternative was measured and rejected: SuperKMeans (arXiv 2603.20009) rotates the data and prunes on a Beta tail, reading 34.1 % of the dimensions on mnist at k=10, but the rotation costs d² per point, which at a CF summary's k is more than the assignment it accelerates. Its own code disables the path below d = 128 or k ≤ 256.

• route_beam — a wider routing descent, off by default. The tree descent commits to one child per level and never backtracks, which is why a quarter to a half of all rows do not reach their nearest microcluster (see Fixed, below). route_beam=b keeps the b nearest nodes at each level instead of one and scans the entries of every leaf the frontier ends on. It reaches every routing entry point — predict, predict_proba, assign_microclusters, outlier_scores, and the sparse CSR path — and nothing else: insertion, the leaf clustering and leaf_refit are untouched, so the tree and the head are identical at every width and only the row → microcluster map moves.

The default is 1, which is the descent this library has always done, so no label changes unless the parameter is set. Measured on ward at max_leaves=4000, median of seeds 0/1/2: the share of re-routed rows converges to the independently measured exact-scan misroute rate on all five datasets (within 0.1–0.8 points), which is what shows the beam finds the exact nearest entry rather than merely a different one. ARI moves on the two datasets whose classes interleave — mnist +0.035 (0.3452 → 0.3801 at b=8) and digits +0.025 at b=2 — and does not move at all on blobs, highdim or covtype, where re-routing 25–44 % of rows lands them in a different microcluster of the same cluster.

The cost inverts the usual trade: blobs/highdim/digits/covtype pay 11–15× on routing at b=16, mnist pays 3.0×, because at d=784 the node features are 24 MB against a 16 MB L3 and the narrow descent is already stalled on memory. End to end that is small: routing is 0.53 s of a 5.82 s mnist fit, so b=8 adds about 0.6 s. The default stays at 1 on the measurement, not on compatibility: b=8 costs 6–10× the routing for zero ARI on three of five datasets, so no single width is right for the table and the choice belongs to the caller. Only heads that assign by microcluster (ward, spectral, leiden, hdbscan) can see the parameter at all, and it is an estimator parameter — the free fit_predict / fit_predict_sparse route with the plain descent. Full sweep in bench/RESULTS.md; CFTree::nearest_entry_beam is public on the Rust side.

• A wheel for free-threaded CPython (3.14t), and a CI gate on what it claims. abi3 cannot express a Py_GIL_DISABLED build — such an interpreter exposes SOABI cpython-314t and no abi3 tag at all — so the single cp311-abi3 wheel matched nothing on 3.14t and pip install betula-cluster there fell back to compiling from source, which needs a Rust toolchain the user did not ask for. maturin drops abi3 for these on its own ("abi3 does not yet support CPython 3.14t … artifacts will be version-specific"), so all that was missing was the build arm: Linux x86-64 + aarch64, macOS arm64 and Windows x64. macOS x86-64 is deliberately absent — the abi3 arm cross-builds it on an arm64 runner precisely because abi3 needs no interpreter at build time, and a version-specific build does.

Measured on cpython-3.14.7+freethreaded before adding any of it: the module imports without re-enabling the GIL, labels are correct, and four Python threads each running a full fit_predict over 6000×8 finish in 0.26 s against 0.22 s for one — 4× the work for 1.2× the time, alternated A-B-A-B. The GIL check is now a CI assertion rather than a note, because a #[pymodule] declaring Py_MOD_GIL_USED turns free threading back off for the whole process on import and does it with a warning rather than an error; a wheel that did that would be worse than no wheel. The Free Threading :: 2 - Beta classifier is deliberately not 3 - Stable: this is one interpreter on one platform, and no test yet drives the extension from several threads at once.

Changed

• The Rust crate's public surface is curated, and the crate now says which modules the contract covers. src/lib.rs declared 23 pub mod, so kernels, linalg and stats — a SIMD reduction, a Cholesky, a lgamma — carried the same semantic-versioning promise as tree and model, and a snapshot of that surface would have been a snapshot of the implementation. Ten are now #[doc(hidden)]: adwin, assign, fidelity, kernels, linalg, mixture, sketch, stats, topology, wasserstein. Nothing is deleted and nothing stops compiling — they stay pub, because this repository's own benchmarks call into kernels and assign — but they are off docs.rs and off the contract, which a crate-level doc section and a new Compatibility section in README.md now state outright. The thirteen that remain documented are the library: tree, feature, distance, bregman, model, clustering, types, sparse, order, coreset, stream, window, validity. What the hidden modules implement stays reachable through the Python package, which is the supported way to use it — mapper() is topology, KllSketch / DdSketch are sketch, the drift distance is wasserstein, tree_report()'s fidelity number is fidelity. Two rustdoc warnings went with the change (40 → 38), and the crate docs gained a compiled minimal-fit example.

The contract is checked by a tool rather than by a hand-kept list: cargo semver-checks --baseline-version 0.8.0 classifies every difference against the published crate, and CONTRIBUTING.md now names it as the gate for a PR that touches a documented module. Run against this tree it reports six breaking changes since 0.8.0, all of them deliberate and all of them already in this section: insert / predict renamed to their try_ forms on CFTree, Model, DenStream, DbStream and WindowStream; the arity of Model::fit, CFTree::build_sharded, clustering::spectral and clustering::kmeans_auto; and this hiding. That is what makes this release a major one.

• save gzip-frames the file, and what load promises across versions is now written down. A model was bare CBOR, and the bulk of one is a leaf centroid matrix that repeats heavily down the leaves — dead dimensions, shared zeros, near-duplicate centroids. Framing it in gzip takes 2.0–4.8× off a real model: digits 0.18 → 0.09 MB, covtype (581 012 × 54) 1.06 → 0.22 MB, mnist (70 000 × 784) 14.35 → 5.29 MB spherical and 22.28 → 9.44 MB diagonal, the largest of them in 0.6 s. Pickle goes through the same bytes, so pickle.dumps and joblib.dump shrink with it, and the container is a plain gzip member rather than a private framing: gzip -dc model.betula still yields the CBOR.

The codec was picked on the artefact rather than on a fragment. zstd -19, brotli -q9 and xz -6 all beat gzip on the same four files — by 12–22 %, for 4–13 s per save against gzip's 0.6 s, which is the wrong trade for a call a user makes interactively — and flate2's default backend is pure Rust, so every wheel target still cross-compiles without a C toolchain. The ratio is a property of the data and is published as one: a standardized Gaussian fixture with no dead dimension compresses 1.1–1.4×, and saying "2–5×" without that row would be quoting the fixtures that flatter the change.

SCHEMA_VERSION stays 2. The CBOR inside is unchanged, so a bump would have refused every existing file to gain nothing; load reads both framings and still rejects any other schema version. That rejection is now a contract rather than an implementation detail — stated in the save/load docstrings, in a new Compatibility section of README.md and in a new Saving and loading a model section of docs/USAGE.md. Models are not portable across major versions, and the alternative — a reader for every past layout — would pin the internals of the CF-tree for the whole 1.x line.

The column-major leaf layout planned alongside this is rejected, and the measurement that motivated it no longer reproduces. On a raw centroid array, column-major moved gzip -6 from 3.97× to 6.15× (mnist at 2000 leaves, 2026-09-08). On the file save actually writes it is worth 1.06–1.17×, because CBOR tags every float with its width and that byte between values dilutes the column similarity the layout exists to exploit. Re-measuring the array itself now reads 2.88× row-major and 3.35× column-major, not 3.97 and 6.15 — the rebuild policy changed which siblings merge, so different centroids end up adjacent. Encoding each per-dimension vector as one CBOR byte string instead of a float array was measured on the same files and is rejected too: 8 %. Neither buys a permanent format complication.

• A rebuild now merges the cheapest sibling pairs, not the closest ones — the leaf summary is 21 % tighter and every published label may move. CFTree::rebuild picked each entry's partner under the routing measure and ranked the pairs by the absorption gap. Under the default D0/D0 geometry both are the unweighted squared centroid distance, so compaction merged the geometrically closest pair whatever it weighed — and the exact price of fusing two cluster features is the Ward cost w_a·w_b/(w_a+w_b)·‖μ_a−μ_b‖², which prices a pair of thousand-point leaves a thousand times above a pair of singletons the same distance apart. Ranking by distance therefore spent the budget on the dense core and left the rest as ones: on digits at 200 leaves 63 % of the leaves held a single point, leaf-mass Gini 0.82. The rank is now the merge cost, whatever distance and absorption are set to; the threshold still grows in the absorption measure, which is the one the gate reads.

Measured against ELKI 0.8.0's BetulaLeafPreClustering at a matched leaf count, the geometry this library defaults to goes from 8 % behind to 15 % ahead: leaf WCSS on digits 7.087e5 → 5.601e5 and on covtype-50k 8.642e4 → 6.794e4, both −21 %, with the singleton share 0.626 → 0.139 and 0.309 → 0.148 and the rebuild count 23 → 8 and 27 → 5. ELKI's own D4/R geometry moves 1.8 % and 5.6 %, which is the control: under D4 the gap already was the Ward cost, and D4/R is the geometry where the two implementations were already level.

Downstream, the ward head gains the most (moons ARI 0.148 → 0.599, varied 0.433 → 0.705, mnist 0.367 → 0.394), gmm gains a little on the real sets, and kmeans loses two real-set cells and gains none (mnist 0.315 → 0.285, covtype 0.064 → 0.049 — both inside the seed spread of the cell). Synthetic k-means and gmm cells move by ≤ 0.001. absorb="chi2" is bit-identical, because its gate binds long before the leaf budget and no rebuild ever runs. Every quality table in bench/RESULTS.md was re-measured for this. (The HDBSCAN contest section was not, and went on quoting the earlier tree; the re-measurement is in [Unreleased].) canonical_order=True is unaffected as a guarantee: the order study rewrote 50 of its 81 cells and all 27 canonical cells still read spread 0.0000 and pairwise ARI 1.0000, which bench/insertion_order.py asserts rather than reports.

Three tuning knobs lose most of what they were worth, because the tree now does their job. balance="auto" fires in none of the 27 cells it was shipped on (the heaviest-leaf share it arms on reads 0.005–0.046 there against up to 0.594 before) and only reaches its trigger two decades of budget lower. leaf_refit=1 now pays in 6 of 12 published column-pairs against 11 of 12, mean effect +0.005 on the arrival order and −0.001 under canonical_order against +0.032 and +0.082 — every leaf_refit=0 cell rose instead. tree_report's size-imbalance diagnosis no longer fires on the fixture it was cut on. All three are documented with the new numbers rather than left quoting the old ones, and none of them changed: the tree they measure did. Full write-up, including the five non-default geometries and the ten-seed spectral control, in research/RESULTS-estep.md.

• The speed half of the benchmark is re-timed on the same tree, and the headline number moves from 8.7× to 10.8×. The rebuild change removes work — rebuilds fall 27 → 5 on covtype-50k and 23 → 8 on digits — and the previous edition published that as an argument for keeping stale timings. It is now a measurement: results_scaling, results_memory, results_real_scale and results_sparse re-run under the page's own contention gate with zero contention sightings, and betula's fit time at a million points is 0.763–0.816× the 2026-09-07 edition across five heads while the scikit-learn column moves 0.952–0.984×, which is machine drift and the control. The effect grows with N: at 10 000 points three of the five heads are slower. 1 M k-means reads 0.216 s, so README.md, docs/MATH.md and paper/paper.md now say 10.8× against KMeans and 37× against Birch; full covtype is 1.006 → 0.746 s at ARI 0.070 → 0.087.

Resident memory did not move at four significant figures, which retires the +7.9 MB covtype regression the 0.8.0 edition recorded: 905.3 → 904.9 MB against an unmoved 931 MB, agreeing with the eight paired runs that already refuted it. The control on the whole re-run is exact — over the 117 rows of results_quality.csv + results_real.csv every metric column reproduced bit for bit and only time_s moved. The scoreboard moves speed 32/1/5 → 30/1/7 and memory 28/6/0 → 30/4/0; both demotions are the same 10 M row counted twice, where betula got faster (8.26 → 7.66 s) and scikit-learn's one-shot path returned to the allocator artefact this page has withdrawn twice (8.96 → 6.43 s, faster at 10 M than at 5 M). They are accepted as measured and named as unreliable in bench/RESULTS.md.

• min_samples now defaults to the leaf mass, not to a point count. HDBSCAN*'s core distance is the radius enclosing min_samples points, and over a leaf summary a single leaf already holds N / max_leaves of them at one coordinate — so any count below that mass is enclosed at radius zero, every core distance collapses and mutual reachability degenerates to plain distance, i.e. single linkage chaining through every overlap. The docs have stated that trap since 0.6.0 and left the caller to avoid it. min_samples=None (the new default, spelled "auto" if preferred) asks instead for the mass of ten average leaves — the conventional ten-point default translated into the currency the head counts in.

On the 100 000-point blobs contest that moves the untuned answer from ARI 0.478 to 0.820 at 2 000 leaves and from 0.678 to 0.896 at 8 000, against fast_hdbscan's 0.910 over raw points. On the six published quality fixtures it costs nothing where the fixed ten already worked (moons 0.9999 → 0.9995, circles 1.0000 → 0.9999, highdim unchanged, aniso 0.568 → 0.565) and gains where it did not (blobs 0.142 → 0.444, varied 0.479 → 0.548). Over a wider grid it wins or ties in 16 of 18 cells, the two losses being a two-blob fixture at a 200-leaf budget where both answers are ARI < 0.08. (Measured before the cost-ranked rebuild above, which shipped in the same release; on 1.0.0 as tagged the rule reads 0.821 and 0.822. The re-measurement, and what else it moved, is in [Unreleased].)

Labels change for method="hdbscan" / "dc-center" / "dc-median" fits that did not name a min_samples; an explicit integer is used exactly as before. get_params() reports None for the automatic state, and auto_min_samples / AUTO_MIN_SAMPLES_LEAVES are public on the Rust side.

• The Hamerly/brute-Lloyd equivalence claim is now true at a cutoff, and tested there. The accelerated k-means is documented as producing output identical to brute Lloyd, and the test that checked it ran only at max_iter = 100, where both have converged. They pair the steps differently — brute is [E, M] × iters and returns the labels of the previous centres, Hamerly is E + [M, E] × iters and returns the current ones — so at max_iter = 1 they disagreed. The reference implementation (test-only) now ends with an assignment step, which is also scikit-learn's contract and the one predict assumes, and the differential test runs at max_iter ∈ {1, 2, 3, 100} comparing labels, centres and inertia, plus asserting that the labels returned are nearest to the centres returned. No shipped path changes: lloyd_hamerly already had the consistent pairing.

• consensus measures one nuisance at a time, and its vote is a bijection — scores move. Two defects, both label- and score-changing. It permuted the insertion order and moved the head's seed on every run while documenting itself as measuring the insertion order; vary now names the choice — "order" (new default: permute, hold the seed), "seed" (one order, move the seed), "both" (the old behaviour). And runs were aligned to the reference by sending each cluster to whichever reference cluster it overlapped most, which is not a bijection: a run that splits a reference cluster maps both halves back onto it, so the split scored as agreement while the merge it forces at fixed k was counted once. The bias has one sign — measured at +0.125 on a split/merge fixture (0.750 against the bijection's 0.625) and +0.03…+0.13 over random disagreements. Alignment is now a maximum-weight bijection (a ~40-line Jonker–Volgenant in NumPy, since the wrapper's only dependency is NumPy; tested against SciPy and against brute-force enumeration at k ≤ 6), and a run with more clusters than the reference keeps its surplus as fresh ids instead of folding them onto a matched cluster. Confidence scores from 0.8.0 and earlier were too high and are not comparable with these.

• Rust API (breaking): the point entry points are checked. CFTree::insert, Model::predict, DenStream / DbStream's insert / predict and WindowStream::insert took a &[R] and trusted its length. The SIMD kernels compare a.len().min(b.len()) coordinates, so a row one column short was clustered on its prefix and a row one column long dropped its tail — both silently, and a debug_assert is not a release check. Each now has a try_* counterpart returning Result<_, ShapeError> (new, in types, with expected / got and a Display), and the unchecked bodies are pub(crate): the Python boundary validates a whole array once and then calls them per row, which is the one caller that should not pay per point. The tree is untouched when a row is rejected, which an integration test pins in both profiles.

For a Rust caller the fix is tree.insert(p) → tree.try_insert(p)?. CFTree::dim() is now public, since a caller that has to match the dimension should be able to read it. The truncation in kernels.rs is no longer load-bearing — every internal caller slices exactly dim — and its comment says so rather than claiming an insert relies on it.

• The streaming constructors bound every parameter they carry, not the product of two of them. DenStream::new tested eps, lambda and beta·mu with is_nan, which an infinity passes, and never tested beta or mu on their own. So eps=inf was accepted and merged the whole stream into one micro-cluster; beta=2.0 set the promotion weight above the macro floor it is defined as a fraction of; and beta=-0.5, mu=-4.0 passed the product test with beta·mu = 2 while leaving mu negative, and mu is the macro-cluster weight floor — every connected component, down to a single stray point, was then reported as a cluster instead of as noise. Each bound is now is_finite and each parameter is checked on its own (eps > 0, lambda > 0, 0 < beta <= 1, mu > 0, then beta·mu > 1); DbStream::new likewise for r, lambda, alpha and min_weight. The rejection reaches Python at the first partial_fit, where the model is built.

• max_iter=0 is refused, and one iteration is the floor everywhere it was not. A zero-iteration EM returns its own initialisation: for the GMM heads a responsibility matrix of zeros and a labelling of all-zero, reported as a successful fit of one cluster. Nine head loops ran 0..max_iter and five of them clamped to one iteration first, so the same request meant different things depending on method. The floor now lives once in model::fit_head, where every parametric head passes, and the three heads dispatched outside it (scale-space, COP-KMeans and the Bregman pair) carry it at their own entry with a comment saying why. On top of that, max_iter=0 raises at every boundary a caller can reach: Betula, fit_predict, fit_predict_sparse, BregmanBetula, KPrototypes and the CLI's --max-iter.

• A head-specific keyword on a head that ignores it now raises. rank=8 with method="ward", min_samples=25 with method="kmeans", resolution=1.5 with anything but Leiden — each configured nothing and said nothing. Nine keywords are now checked against the head that would read them: min_samples / graph_degree (hdbscan, dc-center, dc-median), min_cluster_size (hdbscan), resolution / covariance_weight / tangent_weight / tangent_rank (leiden, leiden-cpm), rank (mppca, mfa) and fuzzifier (fuzzy-cmeans), with the error naming the heads that do read the keyword.

A keyword left at its default is not a request — rank=2 is what every caller who never thought about rank passes — so only a changed value is an error, and the three defaults now live in one place (DEFAULT_RANK and its siblings) that the PyO3 signatures, the serde defaults and the check all read. max_iter and auto_k_max are deliberately not checked: they are general to the parametric heads and their applicability depends on n_clusters as well as on the head.

• A covariance head with feature="spherical" is now refused, not warned about. Spherical keeps one scalar of within-leaf scatter, and cov_dense hands it back as ssd/(w·dim) · I — so a head reading a per-component covariance gets the same isotropic term added to every component, ln|Σ_c| moves with it, and the maximum-posterior argmax follows. On digits at ×2.0 compression method="gmm" reads ARI 0.0097 against 0.4343 on feature="diagonal"; nearest-centre labelling of the same fitted model reads 0.5288, which is how the head rather than the fit was identified as the loser. 0.8.0 warned on the gmm pair only, on the grounds that the other heads were unmeasured. They have now been measured, over five leaf budgets on digits (medians of seeds 0/1/2): gmm-full collapses at two of the five (0.0096 and 0.0115 against 0.6220 and 0.5131) and mfa at ×2.0 (0.0086 against 0.4949), while mppca — whose model is a low-rank subspace plus its own isotropic noise — and the three gmm-toeplitz rungs, which read covariance off the between-leaf structure, are unaffected to four decimals.

So method ∈ {gmm, gmm-full, mfa} with feature="spherical" now raises a ValueError naming the measured cost and the fix, and tune()'s default space stops proposing the pair. An answer that is right at three budgets and near zero at two is worse than no answer, and the caller cannot tell which one they got. bench/leaf_budget.py, bench/size_imbalance.py and bench/insertion_order.py measured their gmm column on that pair; they now pass "diagonal", and those published tables will move when next re-run.

• Rust API (breaking): Model::fit, kmeans_auto and spectral take the restart count. It sits after max_iter, and 0 asks for KMEANS_N_INIT — the value they used before — so the fix at every call site is to pass 0. spectral's private N_INIT and the Python layer's COP_N_INIT are gone, both being the same constant under two names. No fit changes.

• The streaming heads now have a reference, and it says the drift row is a win and the split row is a loss. bench/stream_reference.py runs DenStream / DBSTREAM / CluStream against River's implementations of the same three papers on a 4000-point stream (River is pulled per invocation and is not a dependency). At a shared radius the comparison is a units mismatch — every one of the four density heads reaches ARI 1.000 on the stationary fixture at its own radius, and they simply disagree about what a radius of 1 means. At each head's own best radius, chosen on the stationary prefix: a +50 translation is followed exactly here (1.000) while River's DenStream and CluStream both fall to 0.000, and a cluster split is resolved by River (0.54–0.57) and not by us (0.120 / 0.288), because the radius the settled stream recommends is wider than the halves it has to separate. Learning costs 0.1 µs/point (DenStream) and 1.0–1.7 µs/point (DbStream) against River's 4.7 and 21–23, and CluStream's 150–181. CluStream's pyramidal time frame — retrospective queries over an arbitrary past window — has no counterpart here at all, and bench/RESULTS.md records that as a missing capability rather than omitting the row.

• The ELKI cross-check now covers the tree, and at equal leaf count the default geometry is 8 % behind. cross_check.py always ran a tree-only layer and it was never written up. Read at an equal maxleaves budget it flatters this library — both implementations undershoot the budget and ELKI undershoots harder (926 leaves against 1583 at a budget of 2000 on covtype/D4/R), so its lower WCSS is bought with leaves rather than with a better summary. Binary-searching ELKI's budget until the realised counts match to within 1 % splits the answer by geometry: on D4/R the two trees are level (ratios 0.977 and 0.993), and on D0/D0 — this library's own default — ELKI's leaf partition reaches 7.6–8.5 % lower within-cluster sum of squares. The partitions agree with each other at only ARI 0.39–0.53 even where their WCSS matches, so a matching objective is not a matching tree. research/RESULTS-estep.md records both tables and the timing column that must not be read as a speed comparison (ELKI's seconds are JVM start plus CSV parsing).

• The ELKI cross-check of the GMM head is re-run at five seeds, and it costs the E-step page one cell. research/RESULTS-estep.md claimed the shipped head "leads at the median in all four cells" of an ELKI 0.8.0 comparison; three seeds could not separate a 0.02 ARI gap from its spread. At five, the shipped head's own medians are unchanged (0.5239, 0.5210, 0.0852 — its spread was already resolved) but ELKI's moved, and covtype/D0/D0 is now 0.0632 against BetulaGMMWeighted's 0.0681. That is a tie inside a −0.0168…0.1039 seed range, not a loss, but it is not a lead and the sentence is corrected rather than left standing. The within-cluster sum of squares still favours the shipped head in all four cells, which is the claim worth keeping: WCSS is the objective and ARI is a proxy for it.

The same run makes the two come apart visibly on covtype, where ELKI's k-means scores higher ARI (0.0951 vs 0.0820) while reaching a higher within-cluster sum of squares in every cell — a better-optimised k-means that predicts forest-cover classes worse, because those classes are not k-means-shaped. And one number for the open n_init question: ELKI's 4 restarts buy +0.137 ARI on digits (0.6856 vs 0.5482) and nothing on covtype.

• The uniform 1–12 % slowdown across the scaling and streaming suites is measured, and it is not the insert path. bench/RESULTS.md recorded the shift honestly but could not attribute it, because the compiler moved with the code (rustc 1.98.0 → 1.98.1, kernel 7.1.9 → 7.1.13). Two A/Bs, each varying one thing: the toolchain costs 8.9 % on the kmeans probe, 4.1 % on ward and 7.7 % on streaming, with gmm a wash; the source (0.7.0 → 0.8.0, compiler pinned at 1.98.0, 15 clean samples per cell and 0 dirty of 120) costs 4.1 % on ward and ~2.5 % on streaming, while kmeans and gmm flip sign between cycles and are a wash.

Both probes run bench/_worker.py directly, so the cell measured is the cell published rather than a proxy, and each arm's compiler is read back out of the built .so instead of trusted from the environment variable meant to select it. The second arm nearly measured something else entirely: a fresh uv venv in the baseline worktree resolves the repo's 3.14 pin to uv's managed free-threaded build, not the system CPython the repo's own venv uses, and would have compared two pyo3 configurations rather than two revisions. Both arms now assert the interpreter as well.

What the source change did buy: this edition is 4–12 % lighter in peak RSS on all four probes (88.1 vs 91.4 MB, 88.2 vs 95.5, 52.0 vs 54.8, 52.9 vs 60.1). The kernel could not be varied and so is not excluded, only unnecessary; the reading that is excluded is the one the old note pointed at, since kmeans — the purest insert-bound probe of the four — is 1.5 % faster on the new tree.

Fixed

• A corrupt model file loaded, and the panic arrived three calls later. load checked the schema version and nothing else, so any CBOR document that deserialized became a Betula — including documents describing a tree no insert could have produced. That is the cause: a schema version says which fields exist, not that their values address each other, and a file is input. Of seven structural corruptions built by hand, four panicked: a root index past the arena, a child index past it, a parent index past it and an empty arena all reached index out of bounds inside the next partial_fit, and an entry whose mean was one coordinate wide reached it inside cluster_centers_. A panic crossing the FFI boundary surfaces as pyo3_runtime.PanicException, which is not a ValueError and is not something a caller can be asked to catch.

load now validates every tree it deserializes and raises ValueError: corrupt model: … naming what disagrees. CFTree::validate checks that the arena is non-empty and root addresses it, that every parent and child index is in range and points at the right arena, that no node is reachable twice — which is what rules out a cycle — that every feature is self-consistent and carries the tree's dimension, and that branching, leaf_cap, max_leaves, threshold, huber_k and balance cannot leave an insert unable to terminate. All seven corruptions now fail at the boundary that produced them. The check is a no-op for a tree this process built, which a test pins across fresh, rebuilt, sharded and empty trees and all four feature models — a validator that rejects a real tree would be worse than none.

Two things it deliberately does not do. Unreachable nodes and unreferenced entries are not an error: compaction leaves both behind and they cost only memory. And the decompression bomb the gzip framing invites was checked rather than assumed — decode streams out of the gzip member instead of inflating to a Vec first, and a CBOR length header claiming 2^40 elements is refused without reserving for it. Over 4 000 randomly corrupted files, every failure is a ValueError and none is a panic.

• An f32 tree stopped counting rows at 2^24, and its reported variance drifted for every row after that. A cluster feature accumulates w ← w + wᵢ, and binary32 spaces its values 2 apart from 16 777 216 upward, so a unit-weight row left the weight exactly where it was. Nothing else in the summary stopped with it: the mean went on moving at 2^-24 per row and the scatter went on accumulating against a weight that no longer grew, so the variance drifted with the fraction of rows the weight had dropped — 1.0058 on unit-variance data 200 000 rows past the ceiling, against 1.0000 in f64. microcluster_weights_, cluster_sizes_ and every mass-weighted head read the frozen total.

The weight is now an f64 running total whatever the tree's element type, in all four cluster features, and spherical's scalar scatter with it. That is the split the arithmetic asks for rather than a concession to f32: a weight's precision is set by how many rows a leaf absorbs, which the caller does not choose, while a mean and a per-coordinate scatter are bounded by the data, which the caller picked the dtype for. Both widened fields are O(1) per leaf, so an f32 tree still costs half the memory of an f64 one. What a caller reads back is rounded to the tree's own type on the way out — one part in 16 777 216 of a total that is itself exact. The UserWarning that used to name the approaching ceiling is gone with the defect it warned about.

Nothing already measured moved, and that was checked rather than argued. Across 48 cells (2 fixtures × 2 dtypes × 6 feature/head pairs × 2 leaf budgets) no f64 cell changed at all and no labels changed anywhere, f32 included; the only movement is heaviest_leaf_width in 18 of the 24 f32 cells, by at most 6.5e-7 relative, since the width is derived from the totals that widened. Every published table in bench/RESULTS.md — all f64 — therefore stands unchanged. The on-disk format needed no bump either: CBOR tags a float with its width and serde's f64 visitor accepts an f32 value, so SCHEMA_VERSION stays 2 and a model written by the released 0.6.0 wheel still loads, which a test asserts. Insert-path cost, three sessions per arm alternating builds on 200 000 × 32 at max_leaves = 2000: f32/spherical +1.5 %, f32/full +1.8 %, f32/diagonal and all three f64 cells unchanged.

One residual is left in deliberately. The per-axis and full-matrix scatters of diagonal, full and fd stay in the tree's element type, so past 2^24 rows in a single leaf the scatter is now the term that saturates and the variance comes back 0.89 % low where the defect had it 0.58 % high. Widening those was implemented and measured before being rejected: +11 % on the insert path and an f32 tree three quarters the size of the f64 one instead of half, charged to every f32 user to correct a leaf that has to be past 16.8 M rows before it is wrong at all. an_f32_diagonal_leaf_past_the_ceiling_still_loses_scatter pins the figure and docs/USAGE.md states it.

• leaf_refit handed back leaf statistics for a partition the tree it returned could not produce, and at a large budget that misrouted 15 % of the rows. The pass routes every row through the finished tree and rebuilds each leaf CF from exactly the rows it won — then, if any entry won none, it dropped those and rebuilt the node structure from the survivors. That rebalance replaces the tree that did the routing, so the CFs describe a partition predict no longer returns: a row descends the new structure and is handed an entry accumulated from different rows. Measured on 4 000 8-D points at max_leaves = 160, the routed codebook error read 26.54 against an exact 4.72, 14.8 % of rows on the wrong entry, p99 244 — while the entries themselves were fine, a beam-4 route reading the exact answer with no misroutes at all.

A rebalance is now always followed by another route, so the statistics and the structure describe the same partition. Over 150 cells (5 budgets × 10 blob fixtures × leaf_refit 1/2/3) the routed error rises by more than 5 % in 2 of them against 8, worst 1.33× against 5.25×, and the worst exact nearest-centre ratio falls 1.32 → 1.13; on the cell that blew up every one of the 107 surviving prototypes now wins a row, where the un-refit tree leaves 19 of 144 winning none. bench/results_refit.csv is re-run: the pass pays in 6 of 12 column-pairs with mean +0.006 (arrival) and +0.009 (canonical_order), against +0.005 / −0.001. Labels move wherever leaf_refit > 0, and a pass that drops an entry now costs two routes rather than one — measured at 1.8–2.5× a plain fit for leaf_refit=1 against 1.4–1.8×.

How many discretionary rebalances to take is a cap on work with both extremes measured, recorded on REFIT_REBALANCE_ROUNDS in src/tree.rs: never rebalancing lets a shallow tree drift until its own descent cannot reach its prototypes (39 → 12 entries, exact error tripled, on 1 of 50 fixtures), and rebalancing on every drop cascades to 59 routing passes and prunes 314 entries to 120. One is the best of the five policies measured on every column.

• projection="weighted-nmf" no longer depends on how the threads interleaved. The transpose-product WᵀX is the projection's hot loop and was summed with a rayon fold/reduce, which splits by work-stealing and merges in completion order. Floating-point addition does not associate, so the result was a function of the run, not of the input: measured on a 19 998-row fixture, 20 runs at RAYON_NUM_THREADS=8 produced 20 different components_ bit patterns, while the single-threaded run was stable. The sum is now taken over fixed 64-row chunks and merged in index order, which is one order for one input at any pool size, and a subprocess test pins it. The documented promise that the thread count does not enter the answer was true everywhere else and is now true here; components_ values shift in their last few bits against 0.8.0.

• Every bench/ harness now loads its real-data fixture from one module, and the Zador control stops describing data no cell in its own table had seen. A "20 000-row draw from covtype" was spelled out inline in five separate files, and one-off scripts had written further variants: the quality tables draw with rng(seed).choice and standardize that subsample, while the order and budget studies standardize all 581 012 rows and then take a fixed rng(0) permutation. The same ward head reads 0.0861 under the first rule and 0.1416 under the second, and bench/RESULTS.md printed both under the label covtype-20k. bench/_fixtures.py is now the only definition of any of them and every harness loads through it; the rules are named (sub / pop / raw), tabulated in bench/RESULTS.md, and the three sections that publish a sub fixture under a 20 k label say so in place. All five rewired loaders return bit-identical arrays to the ones they replace across 16 (dataset × seed) cells, so nothing published moved for the refactor.

One number did move, and it was a defect: bench/zador_fit.py's TWO-NN cross-check built its own subsample — unstandardized, rng(0).choice against the sweep's rng(0).permutation, and MNIST taken off the front of the file rather than drawn at all — so it estimated the intrinsic dimension of data the budget sweep never summarised. It now loads through the loader bench/leaf_budget.py itself fits: 4.61 → 5.28 on covtype, 12.91 → 12.94 on digits, 18.65 → 18.50 on MNIST. Re-fitting that table against the results_budget.csv the rebuild change rewrote also retires a published reading: digits ward was the shallowest slope of the four (d_eff 3.67 ± 0.48 against euclidean's 2.85 ± 0.42) and is now the steeper of the pair (3.47 ± 0.39 against 3.73 ± 0.43) — under one standard error either way, so the fit never resolved a routing effect on digits.

• The three studies that swept gmm on the refused feature are re-run, and one published conclusion changes sign. bench/results_budget.csv, bench/results_imbalance.csv and bench/results_order.csv had their gmm columns measured with feature="spherical", which the library now rejects; all three were re-run on 2026-09-09 with the harnesses' current feature="diagonal". Only gmm cells moved — 76 of 92 in the budget study against 0 of the 160 others, which is the control that says the rest of the record is reproducible to the last digit.

The conclusion that changes is canonical_order's cost. It was published as mean +0.0136 over 27 cells, "positive only because of two cells where it rescues a head the arrival order was collapsing" (digits, 360, gmm 0.1738 → 0.5146; mnist-10k, 1000, gmm 0.0551 → 0.2457). Both were the isotropic collapse, so the flag was being credited for rescuing a fit that is no longer offered. On the diagonal feature those cells read 0.4721 → 0.4195 and 0.2223 → 0.2055 — the arrival order ahead in both — and the grid-wide mean is −0.0064, median −0.0052, non-negative in 9 of 27 (was 10). The invariance itself is untouched: 27 of 27 canonical cells still read spread 0.0000 and pairwise ARI 1.0000. docs/USAGE.md and bench/RESULTS.md carry the corrected numbers and say what they replaced.

The budget study's other half is the refusal's confirmation: on the diagonal feature gmm does not collapse anywhere (digits at ×2.0 0.0088 → 0.4403, MNIST at ×5.5 0.0618 → 0.2850), and it becomes the strongest head on digits at ×4.2 and on MNIST below ×20. The spherical numbers are kept, dated, as the evidence for the refusal rather than as current results.

• canonical_order refuses more rows than a u32 rank can address, rather than wrapping. The permutation is a Vec<u32> — four bytes per row instead of eight, which is 17 GB of index at the ceiling — and (0..n as u32) truncates silently past 2³² − 1 rows: the result is a valid-looking order over the wrong rows, with the tail of the matrix never visited. Both entry points now return an error there, canonical_permutation alongside its CSR twin, and it reaches Python as a ValueError. Rust API (breaking): order::canonical_permutation returns Result<Vec<u32>, &'static str>, matching canonical_permutation_csr, and order::MAX_ROWS is public. The arrival-order build is unaffected and has no ceiling.

The row order is now computed once at the Python boundary and handed down, rather than inside the build: balance="auto" builds the tree twice to decide, and was paying for the O(n · dim) projection and the O(n log n) sort on both passes. No label moves — the order is a function of the data, so computing it once and computing it twice give the same ranks.

• Three documentation claims that were no longer true, and a check so the counts cannot rot again. The README and the JOSS paper quoted a 457-case Python suite and 728 Rust tests against an actual 547 and 807, and the README's two sentences disagreed with each other; the counts are now derived by scripts/check_test_counts.py, which collects both suites without running them (cargo test -- --list, pytest --collect-only) and fails on any published number that has drifted. The README's "always faster — the unconditional win … this holds for every method at every size" is now scoped to what was measured — faster and lighter than scikit-learn at every published size and budget — with the two rows we lose named next to it: scikit-learn is faster on raw-TF-IDF 20-newsgroups at a better ARI, and fast-hdbscan reads 0.910 against our 0.478 on 100 k blobs. The paper's "every benchmark figure is the median of three seeds" is true of the quality tables only, and now says so; its mutation-baseline sentence no longer claims an argument for every surviving mutant when what exists is an argument for every recorded one. Betula.save's docstring said bincode; the format has been CBOR since ciborium replaced it.

• The symmetric eigensolver's convergence test is relative to the matrix, not to 1e-15. jacobi_eigen's fixed absolute tolerance answered two questions wrongly at once. Below a ‖A‖_F of about 1e-15 it was satisfied before the first rotation, so the untouched diagonal came back as the spectrum: on a 4×4 fixture scaled by 1e-18, octave-cli's eig reads [-0.4755, 3.3895, 4.0718, 7.2641] after rescaling and the old code read [-0.890, -0.812, -0.440, -0.114] — four wrong values, no error, no warning. At unit scale it was unreachable instead: off(A) sums n(n−1)/2 entries, so every one of them has to fall below 1e-15/n first, and the loop always spent its whole 100-sweep budget in f64 and in f32, where eps is 1.2e-7 and the threshold is meaningless. The test is now off(A) ≤ eps_R · ‖A‖_F (Golub & Van Loan, Alg. 8.4.3), which makes the answer invariant under scaling the input by a positive constant, and a golden test pins five scales from 1e18 to 1e-30 against octave-cli at 13 significant digits. A 32×32 f32 decomposition drops from ~520 µs to ~348 µs on the measurement machine (three alternating repetitions); f64 is unchanged within noise. No label moves: spectral and leiden on digits at seeds 0/1/2 in both precisions return byte-identical labels before and after.

• cholesky_lower rejects a non-finite partial sum. NaN <= 0.0 is false, so a non-finite entry passed the positive-definiteness test, sqrt propagated it, and every logdet and Mahalanobis distance downstream came back NaN — a full-covariance GMM component with a NaN scatter looked like a valid factorisation rather than a rejected one.

• threshold="auto" no longer breaks the canonical_order guarantee. canonical_order=True promises a summary that is a function of the row multiset; the automatic threshold is piloted on a bounded subsample, and that subsample was drawn by row position, so a permutation handed the pilot different rows and it converged to a different threshold — the order dependence was back before the tree was built. Measured on a 12-dimensional 4-blob probe, one permutation moved the pilot from 17.184 to 17.256 and the realised tree from 284 leaves to 279, with the labels no longer identical. The pilot rows are now taken evenly spaced along the canonical order, which is a function of the content, so the threshold, the leaf count and the labels are all invariant again. Only canonical_order=True with threshold="auto" is affected, and its labels change; the arrival-order path keeps the uniform draw and is untouched.

• A loaded model reported projection="none" however it was fitted. Betula.load rebuilds the wrapper with cls(**core.get_params()), and the engine's get_params never reported projection, projection_dim or projection_max_iter — so a model fitted with projection="svd", projection_dim=4, saved and loaded, came back describing itself as ("none", 64). The fitted state survived (components_ and predict were unaffected), but get_params lied and sklearn.clone of a loaded model returned an estimator that does not project at all.

The cause is that the engine stores the projection decomposed — rank, an SVD flag, a KL flag, solver sweeps — because the persisted wire format predates the ProjectionSpec sum type, so reporting it means inverting projection_spec() and nobody did. More generally the engine's get_params is a hand-maintained mirror of the struct with nothing tying the two together, which is how these three and the new route_beam each went missing. A test now asserts the engine reports every parameter the wrapper carries, with memory_budget_mb — resolved into max_leaves before the engine exists — as the one documented exception.

• The greedy descent misroutes a quarter to a half of all rows, measured — and the blast radius is a quarter of what it was thought to be. The premise on record was that every label this library returns goes through CFTree::nearest_entry. It does not: fit_predict ends in route_data, which consults the tree only when no PointRule is set, and finalize sets one for every head with a centroid assignment rule. On digits, kmeans labels agree with an exact argmin over the k cluster centres 100.00 % of the time and gmm uses the mixture posterior — neither touches the tree at label time. Only ward, spectral, leiden and hdbscan assign by microcluster.

For those, against an exact BLAS scan at max_leaves=4000: 24.6 % (blobs) to 44.4 % (covtype) of rows do not get their nearest microcluster, and a microcluster's own centre fed back through the tree returns that microcluster only 48–71 % of the time. It costs almost nothing on four of the five datasets — 0.00–0.23 % of labels change, ARI between the two routings ≥ 0.994, and against ground truth greedy is a wash or marginally ahead — because the median misroute lands 0.05 to 1.15 centre spacings away, i.e. in a neighbour that belongs to the same cluster.

mnist is the exception: 9.5 % of labels change and exact routing is worth +0.037 ARI (0.3533 → 0.3902, a 10 % relative gain). At d = 784 the excess is the smallest in the table (0.06 spacings) and still crosses cluster boundaries — concentration of measure makes the wrong answer cheap in distance and expensive in labels at once. No unconditional fix is taken: an exact scan costs 0.7× the fit on mnist but 8.9× on covtype (and covtype's scan is the cheaper of the two in FLOPs, so the cost is n × m memory traffic, worst exactly where betula's advantage lives). An exact scan below an m₀, or a beam at descent, both change mnist labels and so are type-1 doors — the beam is the one taken, as the opt-in route_beam above. Full table in bench/RESULTS.md.

• A published claim retracted: the +7.9 MB peak-RSS move on full covtype is not reproducible. bench/RESULTS.md called it "a real, reproducible move in the wrong direction" on the strength of one sample against one sample — two paragraphs after warning that a single sample cannot settle this cell. Eight paired runs of bench/_worker.py real_fit betula-kmeans covtype, alternating this tree against 28d4e9f, put the medians at 884.7 MB against 882.6 MB — 2.1 MB apart, inside within-arm spreads of 9.6 and 15.3 MB. The same eight pairs make the time gain larger than published, 1.005 s against 1.158 s (13 %) where the table reads 1.06 → 1.01 s. The scoreboard cell is a tie either way; what changes is that there is no unattributed memory regression to chase.

• The weekly mutation run was sharded 96 ways on a mutant count that has since grown 2.1×, and more than half of its shards are being cancelled at the timeout. Found by scripts/check_mutants_baseline.py, which prints the live count as a side effect: cargo mutants --list reads 8 927 against the 4 312 the matrix was sized on. Run 34108190109 was in flight while this was being fixed and settled it — 28 of its first 53 completed shards were cancelled at exactly the 150-minute cap. A cancelled shard still uploads its partial mutants.out, so the collation job counts it as present and the week's survivor total simply comes out short, with nothing red to say so.

The 23 shards of that run that did finish give a cost model, taken from the per-mutant phase timings in their outcomes.json rather than from a per-mutant average — averaging is what hid this, because the distribution spans two orders of magnitude: 15 s unviable, 70 s caught (p90 447, max 1082), 660 s missed, 988 s timed out, plus a 236 s baseline per shard. -j4 measures 0.94 parallel efficiency, so job_minutes ≈ worker_seconds / 4 / 0.94 / 60. Shard cost is set by the mix, not the count: those 23 shards span 7 024 to 34 971 worker-seconds at an identical 93 mutants each, because --shard k/n slices the deterministic list contiguously and so effectively by file.

That rules out raising the cap alone: 93 mutants that all time out is 467 minutes, past GitHub's 360-minute ceiling for a hosted job, so no cap makes a 96-shard matrix safe. The matrix is 192 shards with a 300-minute cap, which puts that same worst case at 236 minutes and the worst observed mix at ~78. Wall clock does not pay for the split — the total work is unchanged and only the per-shard baseline and job setup are duplicated, ~25 minutes against a measured total of ~610 worker-hours. The shard count now lives in one env: SHARDS that both --shard k/n and the collation tally read, so the two can no longer drift apart silently; a mismatch with the literal matrix surfaces as shards collated: 96 / 192 and fails.

The v0.8.0 changelog reached the "cannot finish as configured" conclusion by wrong arithmetic (a stale 101 worker-seconds per mutant, and -j4 forgotten); that sentence is corrected in place rather than deleted.

0.8.0 — 2026-09-07

Fixed

• A column-major (Fortran-order) input array was read as its transpose — silently, with no error, producing a clustering of scrambled data. Every array-taking entry point was affected: fit, fit_predict, predict, partial_fit, and the mixture_w2 parameter readers. It is not an exotic input — sklearn.datasets.fetch_openml returns one, as do np.asfortranarray, several pandas .values paths, and interop with R or Julia.

The cause is one word in the zero-copy ingest (06ad6b5): PyReadonlyArray::as_slice succeeds on any array that is contiguous, and a column-major array is contiguous. Its buffer is then the transpose of what a row-major read expects, so the borrow handed the tree the right numbers in the wrong shape. Nothing downstream can notice: the values are finite, the moments are well-formed, the tree builds, and the labels come back. On MNIST (fetch_openml, 10 000 × 784) it read as one leaf holding 2381 "identical" rows, a summary no data row was within 2273 of, and an ARI against the digit labels of −0.0001 where scikit-learn's k-means reads 0.34. After the fix the same call reads 0.344 at max_leaves = n, with a lower k-means objective than scikit-learn's (2.5238e10 against 2.5284e10).

The borrow is now gated on ndarray's standard layout — the row-major condition the module actually indexes by — rather than on contiguity, and a column-major array falls through to the existing copy path, which iterates in logical order and was always correct. mixture_means and mixture_covs had the same defect and now read logically too. The regression test asserts label equality between the two layouts rather than similarity, and fails on revert.

No published number moves: every dataset in bench/ arrives C-contiguous (fancy indexing and StandardScaler both produce one), which is also why this survived so long — the benchmark record could not see it.

• Two reference values in the Watson head's GOLDEN table were computed in double precision and printed as if they had 40 digits. No shipped behaviour changes; the table is test input, and the test's tolerance (1e-11 relative) is four orders of magnitude looser than the error, so it passed before and passes now. The cause: local/scratch/watson_kummer.mac passed the float literals ±0.25 and 1e4 where exact ±1/4 and 10000 belong. A float argument makes Maxima's hypergeometric evaluate in float mode and return a double, which bfloat then pads out to fpprec: 40 — 15 real digits followed by 25 of noise. Rounded to f64 the d = 3, κ = +1/4 row was 5 ulp wrong (8.615403391270196e-2 against 8.615403391270206e-2) and κ = −1/4 1.6 ulp; both are corrected. The same float literal made section 5 — the section that verifies the bounds g(0) = 1/d, g(−∞) = 0, g(+∞) = 1 that let the concentration solver pick its branch — abort on a floating-point overflow, so that check had never actually run; with exact arguments it runs and confirms all three to 40 digits.
• Six reference values in stats.rs's CDF_TABLE were evaluated at the decimal x rather than at the f64 the test passes to reg_lower_gamma. Also test input only, and also below the tolerance that reads it (5e-13 against the implementation's own worst residual of 1.1e-13), so nothing was ever red. P(a, x) has logarithmic sensitivity ≈ a in the left tail, so the ~1e-17 gap between the decimal literal and the binary64 nearest it is amplified by a: the six affected rows are exactly the two smallest-x entries at a = 12, 50, 200, and the error runs from 2.3 ulp up to 61 ulp at a = 200, x = 10.05. The predicted amplification (60.5 ulp) and the measured one (61.3) agree, and all six committed values reproduce exactly when the reference is recomputed at the decimal — which is what identified the cause.

• A generator that silently drops to double precision leaves the provenance comment true and the numbers wrong, so provenance is no longer accepted as evidence. All three reference tables in the crate are now recomputed from an engine other than the one that produced them (local/scratch/recheck_goldens.py, mpmath at 60 dps against Maxima and against mpmath's own earlier output): watson.rs agrees to 1.19 ulp over 48 values, vmf.rs to 0.48 over 11, and stats.rs's 234-row quantile table to 0.50. Relatedly, maxima -b exits 0 even when the batch aborts mid-file, which is why the Watson overflow went unseen; local/scratch/run_maxima.sh now runs every .mac under -Q, scans for the abort signature and requires a completion sentinel, and all 13 scripts pass it.

• The release workflow would have failed after publishing. release.yml pipes this file's section for the tag straight into the GitHub Release body, and GitHub caps that body at 125 000 characters. 0.6.0's section is 16 649 bytes and published fine; 0.8.0's is 172 863, so the github-release job would have returned 422 — after the publish and crates-io jobs had already put the wheels on PyPI and the crate on crates.io, leaving a published version with no Release and therefore no Zenodo archive, and no way to re-run it without a new version number. The step now truncates on a line boundary at 110 000 bytes and appends a link to the full file; under the cap it is unchanged. Verified by running the step's script against this tree: it emits 108 517 characters out of 170 426.

Changed

• CFTree::build_parallel is now CFTree::build_sharded, is compiled without the parallel feature, and takes a trailing order: Option<&[u32]> — a breaking change for Rust crate users. Python callers are unaffected; the migration is a rename plus one argument, and None reproduces the previous behaviour exactly.

The order parameter exists because sharding had to split ranks of the canonical order rather than rows: sharding rows would have kept the insertion-order dependence while looking like it removed it, which is a failure mode no test would have caught by accident, so there is one that fails on exactly that mistake.

The rename and the un-gating go together, and both come from the same rule — the thread count and the build configuration are not allowed into the answer, only the shard count is. The old name said "parallel" for something whose observable effect is a partition, and the old #[cfg] meant a crate built without parallel silently produced a different summary for the same n_jobs. rayon now supplies the schedule and nothing else: the partition, the per-shard budget and the merge order are identical either way, and the tests that cover them run in both configurations rather than only under --features parallel.

• fit_predict's n_jobs is renamed n_shards, now int | None defaulting to None — a breaking change for Python callers who passed it. There is no deprecated alias: the migration is a rename, and n_shards=1 or the default reproduces the previous default exactly. consensus(n_jobs=…) is untouched — there it really is a worker count over independent runs, and the result is invariant to it.

It was never a worker count here. It was the shard count, and shards are the partition: two counts hold different point sets and build different sub-summaries, which no merge order repairs. Measured over the 27 published cells at 1 / 2 / 4 / 8 shards, labels at 8 agreed with 1 at pairwise ARI 0.46 on average, 0.098 at worst wherever compression was real — as wide as the row-order gap canonical_order exists to close. scikit-learn documents n_jobs as a worker count and no more, so the old name broke a convention rather than a written contract; but the convention is real — no scikit-learn estimator moves with it — and it is what the name buys. Parallelism comes from RAYON_NUM_THREADS, as it already does for the kernels and the Phase-3 heads, and does not enter the answer.

• n_shards together with canonical_order=True raises ValueError rather than being ignored. With the flag on the shard count is clamp(n / 25000, 1, 64) — a function of the data — because a guarantee that survives a reshuffle but not a re-tuned shard count is not a guarantee. Silently discarding an argument it will not honour is the shape that let the old name go unexamined for two releases, so the two are rejected as a pair instead. The constant is the measured knee rather than a round number: at n = 200k, max_leaves = 2000 a kmeans fit reads 1.00 / 1.90 / 2.66 / 3.42 / 3.71 / 3.91× for 1 / 2 / 4 / 8 / 16 / 32 shards, so 25 000 rows per shard buys 87 % of the available speed-up and leaves the partition four times coarser than the alternative. Sharding's quality cost was measured before the constant was picked and is a wash (8 shards minus 1, mean +0.004 ARI). Inputs under 25 000 rows return one shard and keep the plain sequential build, so every published cell is unchanged. With the flag off nothing changes at all.
• The canonical order's "3 → 30 rebuilds, insert 1.56×" was a tail draw published as the cost of high dimension — one claim corrected, no code change, and a magic constant deliberately left alone. order::PROJECTION_SEED seeds the eight projections the insertion key is built from, and the worst cell on the cost page made it look badly chosen: a seven-candidate screen over six shapes found four constants reading 3 → 3..5 at 50k × 784, max_leaves = 8000 against the shipped constant's 3 → 30, and 138–152 total rebuilds against 172.

It is not badly chosen, and the screen could not have shown that either way. Quality does not move with the constant (18 cells × 7 constants × 3 head seeds; Friedman p = 0.53, largest gap between any two constants' mean delta 0.028 against the arrival order's own median spread of 0.077), so the choice was a cost choice — and the cost margin vanished on shapes the candidates had not been screened on (107 / 109 / 99 total rebuilds for shipped / winner / runner-up, with the roles inverting outright on blobs 60k × 384: shipped 3, winner 15). Twenty-four draws per shape give the reason: the rebuild count is a right-tailed lottery over the projection draw, reading min 3, median 4, max 31 at the very shape that started this, with the shipped constant at 4. One draw in roughly twenty-four is pathological at a given shape, and which one changes with the shape.

So the constant stays, and stays untuned: a screen of a few shapes will always produce a winner that does not generalise. Selecting the draw from the data at fit time would have preserved the order-invariance guarantee and is rejected on measurement instead — no cheap statistic of the codes predicts the rebuild count (Spearman |ρ| ≤ 0.48 over those draws, mostly under 0.27, sign-flipping between cases), so there is nothing to select on short of building a tree per candidate. Recorded in docs/adr/005-canonical-order-projection-seed.md. - bench/RESULTS.md's reproduce recipe pinned OpenMP, which hands betula a free 2.2× — recipe corrected, no published number moves. The page told the reader to export OMP/OPENBLAS/MKL/NUMEXPR_NUM_THREADS=1. Following it does not reproduce the page, because the two libraries do not share a thread knob: betula's parallelism is rayon's, which OMP_NUM_THREADS does not reach, while scikit-learn's k-means Lloyd loop is OpenMP. Pinning OpenMP gags one side and leaves the other on every core.

Every table was in fact taken at OMP_NUM_THREADS=8, inherited from the shell — which makes comprehensive.py's os.environ.setdefault a no-op for that variable and leaves it pinning only the harness's own BLAS. Instruction and numbers were never the same run. Re-running the whole sweep with that single variable changed, scikit-learn alone slows 2.2–2.7× wherever d is large enough for the Lloyd loop to matter — 1.03 → 2.80 s on the memory suite at n = 500 000, d = 20, 4.92 → 10.86 s on full covtype, 1.03 → 1.80 s on 20-newsgroups — and does not move at all on results_scaling.csv, whose blobs are two-dimensional. That dimension dependence is what identifies the cause: a dependency bump or machine drift would not spare d = 2. Two alternative explanations were tested and rejected first — scikit-learn 1.7.2 and 1.8.0 are slower than the installed 1.9.0 at equal threads (22 s against 13 s on covtype), and every affected cell's ARI is unchanged to three decimals, so it is not a different computation. - Every table in bench/RESULTS.md re-measured 2026-09-07, and it is the first edition taken under a contention gate. This machine is shared with other long-running jobs, and the two libraries do not lose the same fraction to a competing process — scikit-learn's OpenMP loops lose more than betula's rayon does, so contention moves a ratio while leaving every label and every structural count bit-identical. Repetition cannot detect that. Each of the three measured steps waited for four consecutive clean 15-second samples, was sampled every 10 seconds while running, and was discarded and re-run if any foreign process crossed 80 % of a core; the first took three attempts.

The quality column reproduced and the speed column did not, which is what makes the speed changes attributable. 108 of the 117 rows in results_quality.csv + results_real.csv are byte-identical in every metric column to the 2026-08-24/25 edition across 180 commits — every scikit-learn row, and the whole of results_real_hires.csv, results_real_normalize.csv and results_refit.csv. The nine that moved are all betula-spectral, from the rewrite already described above.

What moved on speed, and why: betula-hdbscan at 1 M 0.801 → 0.517 s and, where those heads dominate, far more — on mnist 20 000 × 784 hdbscan 15.61 → 1.75 s and ward 17.67 → 3.62 s, which is 06ad6b5 putting the O(m²) heads on rayon. Streaming peak RSS fell 60.2 → 52.9 MB, uniformly at all five sizes, because the suite streams 50 000 × 20 f64 chunks — 8.0 MB — and the same commit stopped the ingest duplicating each one; the saving is one chunk and not a function of N. Headline figures move with it: 1 M k-means 0.26 → 0.28 s, 9× → 8.7× against scikit-learn's KMeans, 30× → 29× against Birch, 82× → 94× on 10 M streaming memory. README.md, docs/index.md, docs/MATH.md and paper/paper.md are updated to match.

Three things are recorded as losses rather than smoothed. (1) Every betula row other than hdbscan on the d = 2 scaling fixture got 1–12 % slower, and the streaming time column by the same 6–9 % at every size — a uniform shift across two independent suites is not a timing sample, it is either the insert path or the toolchain (rustc 1.98.0 → 1.98.1, kernel 7.1.9 → 7.1.13), and this run cannot tell which. (2) betula's peak RSS on full covtype rose 905.3 → 913.3 MB against an unmoved scikit-learn, which demotes two scoreboard cells to ties. (3) Every time in the 20-newsgroups table fell, some by a factor of two, with every ARI byte-identical. Three further gated repetitions with the six arms interleaved put five of the six published values above the whole four-sample gated range, which is what a slower machine looks like — but they also put the betula-sparse / sklearn-kmeans ratio at 2.18–3.22 with the published 2.23 inside it, so that ratio is not resolved by four samples and the page now says so rather than quoting it. Within a gated session every arm slows 4–21 % from the first repetition to the third, which is most of that spread. The one pair the samples do resolve is betula-nmf against sklearn-nmf: a tie inside 6 % in all four, against a published 0.79. The board therefore moves 33/0/5 → 32/1/5 on speed and 30/4/0 → 28/6/0 on memory, accepted with scoreboard.py --update in the same commit.

One published cell was a transcription error, not a measurement. betula-hdbscan on MNIST was printed in the real-datasets table as 0.000; its own results_real.csv has read 0.117 in both editions. The table is corrected; nothing was re-run to correct it. - The scoreboard ratchet identified a cell by who won it, so a champion changing read as a result vanishing. bench/scoreboard.py --check is now a CI job, and it could not have been one before: on a clean tree it failed, reporting two results_sparse cells as VANISHED. Neither had moved — betula-svd had overtaken betula-sparse as our fastest sparse head, which renamed the key …/betula-sparse-vs-sklearn-svd out of existence while the verdict stayed loss. The cause is in the producer, not the check: vs-best and vs-external pick a champion per side, and folding the winners' names into the cell's identity made an improvement indistinguishable from a missing comparison. A cell is now identified by what it compares — axis, table, pairing, dataset slice — and the champions are printed on the rendered line. scoreboard.json is migrated in place; only the 44 vs-best / vs-external keys change, and no verdict does. - BIRCH's size-imbalance bug is D0's, not the CF-tree's — two published claims corrected, no code changed. bench/size_imbalance.py reproduced scikit-learn #22854 at absorb="euclidean" only, and concluded the mis-allocation was "a property of the shared design — one global absorption radius". The sweep now runs all eight absorption criteria, and they differ sharply: on structured (one core of 20 000 against six minorities of 30, medians of seeds 0/1/2) only ward (D4) and chi2 hold ARI 1.0000 at every leaf budget; manhattan, diameter and radius recover at 4 000 leaves and not below; euclidean and average never recover — sixteen times the budget moves neither off 0.4174. ward is immune at 250 leaves with the same realised leaf count as the default (247 against 241), so it is not a larger budget, it is the same budget spent by a criterion that does not measure a radius; chi2 gets there with half the leaves. That makes absorb= a second independent fix alongside balance=, and the cheaper of the two — a bigger budget is not a third, since it rescues D1/D3/R and leaves the default exactly where it was. docs/USAGE.md separately claimed the variance-minimising criteria "inherit BIRCH's size-imbalance bug" — D4 is the variance-minimising criterion and it is the best cell in the table; that paragraph now names which criteria are actually affected. bench/results_imbalance.csv grows the absorb rows so both corrections are reproducible from the shipped harness.

• The spectral head clusters every leaf instead of 256 k-means landmarks, and it is 2–12× faster for it. Above 256 microclusters method="spectral" used to reduce them to 256 weighted landmarks, cluster those, and let each leaf inherit its landmark's label — a second lossy summarisation stacked on the one the tree had already done, forced by the O(M³) cyclic-Jacobi eigensolver. The top-k eigenvectors of P = D^{-1/2}AD^{-1/2} now come from Chebyshev-filtered subspace iteration: the iteration only ever needs the product P x, which is O(nnz) on the sparse graph; a degree-12 Chebyshev polynomial evaluated by the three-term recurrence multiplies the wanted end of the spectrum by cosh(12 · arccosh σ) while staying bounded by 1 on the damped interval, and the Rayleigh–Ritz step runs on a (k+4)×(k+4) matrix. Nothing M×M is formed. Past 2048 nodes the affinity graph comes from the bounded-degree beam-search index as well (the same graph the HDBSCAN head runs on), since the O(M²) distance matrix is the next thing to become unaffordable — 134 MB at 4096 nodes, 3.2 GB at 20 000. A/B on identical trees, both sides rebuilt, median of seeds 0/1/2: two-moons and two-circles at N = 20 000 hold ARI 1.000 at every budget from 500 to 5000 leaves while the fit drops from 0.34 s to 0.02 s at 500; digits-PCA20 goes 0.660 → 0.779 at 500 leaves and 0.786 → 0.801 at 1000, and loses one budget, 0.766 → 0.735 at one leaf per point; covtype stays within ±0.01 of zero on both paths, which is its documented failure case. At or below 256 microclusters the code path is unchanged and the labels are identical. Labels change above 256 microclusters. The eigenvectors are checked against the exact solver by subspace angle under Davis–Kahan's sin θ_max ≤ ‖R‖_F / gap, not by ARI — on three fixtures whose eigengaps span an order of magnitude the measured angle tracks the gap from 1.5·10⁻⁸ to 1.5·10⁻⁵, and a fixed tolerance would pass a broken solver on one and fail a correct one on another. The head's remaining ceiling is the graph's, and the control says so: at 20 000 nodes it scores ≈ 0.60, and forcing the exact O(M²) affinity there returns the same answer for 17.2 s instead of 1.7 — what breaks is KNN = 10 falling to 1.0 · log n, the García Trillos–Slepčev connectivity threshold. Tables in bench/RESULTS.md.

On MNIST the change is a budget shift, and the A/B above missed it — this is the loss. That A/B covered two-moons, two-circles, digits-PCA20 and covtype, none of them above 64 dimensions after reduction. At 784 the benchmark's fixed max_leaves = 4000 reads ARI 0.203 → 0.101, with the seed ranges disjoint (published [0.183, 0.203, 0.248] against [0.079, 0.101, 0.122] over seeds 0/1/2) — and out of 117 quality cells, the nine that moved at all are all this head's.

It is neither the new solver nor the approximate graph. Sweeping the budget through both boundaries at a fixed subsample, quality rises monotonically — 256 → 0.0005, 512 → 0.007, 1024 → 0.025, 2048 → 0.050, 4000 → 0.082, 8000 → 0.203 — and the ≤ 256 row, whose code path is unchanged from 0.7.0, is the worst of the six. What the landmark path bought was smoothing: k-means down to 256 weighted landmarks conditions the affinity graph better than 4 000 raw 784-dimensional centroids do, so the rewrite needs roughly twice the microclusters to build the same embedding.

It pays for them many times over, and then passes the old head. At equal quality the new head reads 0.2029 in 6.21 s against the landmark path's 0.2033 in 43.87 s on the same machine and harness one commit earlier — a 7.1× speed-up, measured on a contended machine and therefore a lower bound. Given the resolution outright it does better than the landmark path ever did: the curve continues 12 000 → 0.235 and peaks at 16 000 → 0.328, 61 % above the published 0.203, before turning down at one leaf per point (20 000 → 0.299) exactly as the documented graph ceiling predicts. Sixteen thousand is what bench/RESULTS.md's hires table already gives every head.

The fixed-budget row is published as the loss it is. The correction that matters to users is in docs/USAGE.md, which advised spending "a few hundred to a few thousand leaves" on this head — advice derived from two-dimensional non-convex fixtures and digits-PCA20, and actively wrong at 784 dimensions, where it costs a factor of four in ARI.

• clustering::SPECTRAL_MAX_NODES is now SPECTRAL_EXACT_NODES, and a second constant SPECTRAL_DENSE_GRAPH_MAX joins it. Rust API only; the Python surface is unaffected. The old name described a landmark cap that no longer exists — the value 256 is now the size at which the exact eigensolver hands over, not the size the node set is truncated to, and the two boundaries the head actually has are worth naming separately because they exist for different reasons (O(M³) solve, then O(M²) memory).

• linalg::orthonormalize_rows is public, and mppca / mfa / spectral share it. Gram-Schmidt with one re-orthogonalisation pass existed twice already, byte-identical; the Chebyshev solver needs a third caller, and three copies of a numerical primitive is how they drift apart.

• clustering::xmeans is now clustering::kmeans_auto; the name xmeans belongs to the new splitter. Rust API only — the Python method strings are unaffected. The function was never x-means: it fits a full k-means at every k ∈ [k_min, k_max] and keeps the best BIC, which is a BIC-selected sweep. X-means is defined by its per-centre split test, and now that the crate ships that too, two pub items could not both carry the name. betula_cluster::clustering::xmeans was published in 0.6.0, so this is a breaking rename for Rust callers; method="kmeans" with n_clusters=0 still routes to the sweep and its labels are unchanged by the rename.

• method="scale-space" sweeps the bandwidth in two passes, and stops answering k = 1. The head ranks flat runs of the mode-count-versus-log h curve, and the swept range runs to half the data diameter — far above the scale at which the last two clusters merge. Eleven of fifteen grid points therefore sat in the fully-merged tail, and the longest flat run of a decreasing curve is always its tail: on digits the head returned a single cluster in 42 of 52 (PCA dimension × leaf budget) cells. The sweep now stops at the first single-mode bandwidth, leaving that run one point long, then narrows onto the merge cascade and re-sweeps it once — guarded on having seen at least two multi-mode scales, because refining below the merge scale manufactures structure out of a single blob. Same 52 cells: 46 wins, 3 losses, 3 ties, mean ARI 0.3845 against 0.0710, and 3 cells at k = 1 instead of 42, for ×2.45 wall clock (one cell runs ×0.24, four times faster, because truncation stops it after two bandwidths). covtype remains the head's worst dataset — every cell there scores a negative ARI, which the selector cannot fix. Labels change on almost every input. Tables in bench/RESULTS.md.

• method="scale-space" mollifies each leaf with its own scatter, and the head's operating envelope is now published. The head advertises the density modes of the data but was building the kernel density of the leaf centroids, one point per leaf. A leaf is a cloud, and convolving that cloud with the kernel is what the data's density is: N(μ_j, Σ_j) * N(0, h²I) = N(μ_j, Σ_j + h²I). Every leaf now carries its own width s²_j = h² + σ²_j with σ²_j = S_j/(n_j·d) from the summary, its own s_j^(−d) amplitude — without which a fat leaf peaks as high as a tight one — and the variable-width mean-shift fixed point. At zero leaf scatter this is bit-identical to the previous path, which is a test rather than a claim. Labels can change where leaf radii vary.

Measuring it turned up something larger. Over 52 (PCA dimension × leaf budget) cells on digits the correction is 7 wins, 7 losses, 38 ties — and 38 of the ties are k = 1 on both sides; on covtype-20k every cell either returns one cluster or returns modes at negative ARI. The best cell in the grid (ARI 0.59) sits one budget away from a k = 1 cell. The fragile part is the plateau selector, which counts grid points on the mode-count-versus-log h curve rather than measuring a persistence, so it moves whenever the sweep moves. The mollification ships because it is the model the head documents and costs nothing, not because it measures better. Envelope and tables in bench/RESULTS.md.

Performance

Three ceilings the measurements themselves named, each attacked with the label digest as the gate: a fit that is 90–99 % insertion, Phase-3 heads that were all single-threaded O(m²), and an ingest that duplicated the input. Every number below is a median of three alternating (A-B-A-B) repetitions on one machine (Zen 3, 8C/16T, AVX2, no AVX-512), OMP_NUM_THREADS=1, with the two arms differing only in the code under test.

• The input is no longer copied when it does not have to be. fit / fit_predict / partial_fit / predict read a C-contiguous array of the tree's own dtype in place; a copy is made only for a strided view, the other float width, or a head that prepares rows (normalize=True, the directional heads, method="hyperbolic") and so must not write into the caller's array. On 200 000 × 784 f64 — a 1 196 MB input — peak resident memory falls from 2 485 MB to 1 357 MB, i.e. the duplicate is gone rather than reduced; on 1 000 000 × 20, from 350 MB to 197 MB. Label digests are unchanged on every shape measured.
• Ward and HDBSCAN* use every core. Both heads were exactly O(m²) on one thread, and at large max_leaves that made the head, not the tree, what a fit waits on: at 7 250 leaves of d = 784, Ward took 67 s against a 6.8 s build. Ward's nearest-neighbour scan, HDBSCAN*'s core distances and Prim's relaxation now split across threads. At 3 631 leaves of d = 784, median of 3: Ward 9.76 s → 2.74 s (3.6×) and HDBSCAN* 9.09 s → 1.93 s (4.7×) on 8 threads. Every reduction is exact — an index-ordered collect, disjoint per-index writes, or a min whose tie-break is a total order on (value, index) — so the labels are bit-identical at 1, 8 and 16 threads and identical to the serial release. The argmin inside Prim stays serial deliberately: it is O(m) scalar compares against O(m · dim) distance evaluations, about 20 ms of a 46 s run. Ward's scan also walks a dense list of surviving clusters instead of a mask over all m; the first shape of the parallel scan forked on every scan whatever its length and returned only 1.21×.
• aarch64 gets the SIMD kernels it never had. The packed sq_euclidean / dot / manhattan reductions existed for AVX2 only, so the macOS-arm64 and linux-aarch64 wheels ran the scalar fold. There is now a NEON path, and it needs no runtime detection — ASIMD is mandatory on AArch64 — so unlike the x86 path it can inline into the caller. Verified under qemu-aarch64 against the scalar fold at every width and tail length; the length threshold is inherited from the measured x86 crossover and is not measured on ARM hardware, which is stated where the constant is.
• A profile-guided build of the linux x86_64 wheel, measured but not yet on the release path. bench/pgo_train.py is the training workload and .github/workflows/pgo-wheel.yml the two-pass container build, dispatch-only for now. Insertion at 1 000 000 × 20 runs 0.948 s → 0.841 s (1.13×) with the profile, with no overlap between the arms across three repetitions; at 200 000 × 784 it is a wash (2.643 s → 2.747 s, ranges overlapping) — the win is where the time is branches and pointer chasing rather than kernel arithmetic.
• benches/, with no benchmark framework. cargo bench --bench kernels sweeps the three kernels over d ∈ {8 … 1024} in both widths, and --bench tree_insert times the insert path through the Rust API, so a layout change can be measured without a maturin develop in between. harness = false and std::time throughout: a dependency in [dev-dependencies] still has to be justified, and the numbers that matter here (rows/s, GB/s at a given d) are not the ones a framework reports. What the first run says is the shape of ceiling 1 — the kernels reach ~70 GB/s at d = 64…128, where a pair fits in cache, and fall to 17–20 GB/s at d = 784.

Added

• canonical_order= — the CF-tree summarises a multiset instead of a sequence, and the row order stops mattering entirely. This is the defect the leaf_refit entry below could not close, and it is structural to every BIRCH-family tree rather than a bug in this one: the prototype set is a greedy threshold-net built in arrival order, so which regions get a prototype — and how many — is a function of the sequence. leaf_refit moves each prototype onto the centroid of its cell but never changes which cells exist, which is exactly why it improved quality and left the order gap intact.

With canonical_order=True the rows are sorted by a key derived from the data before the first insert, and the guarantee is exact: two permutations of the same rows produce label arrays that are equal element for element. Pairwise ARI between orders goes from 0.22–0.68 to 1.0000. The key is a Morton code over eight random projections drawn from a fixed constant — deliberately not from seed, so varying seed re-runs the head over one fixed summary — with ties broken on the full row, because a key alone is not a total order and colliding rows would otherwise keep the caller's order and void the whole guarantee. New betula_cluster::order module; the sharded path is invariant too, because sharding splits ranks of the canonical order rather than rows — see the CFTree::build_sharded note under Changed.

It buys reproducibility, not accuracy, and the measurements are what say so. In the published sweep (bench/insertion_order.py, 27 cells, eight permutations each) the change against the arrival order's median draw is mean +0.0136, median −0.0017, non-negative in 10 of 27, and the canonical value falls inside the order arm's own [min, max] in 21 of 27 — usually indistinguishable from a draw you would have got anyway. It fixes which draw you get rather than moving the distribution; the positive mean is two cells where arrival order was collapsing a gmm head (digits@360 0.1738 → 0.5146, MNIST@1000 0.0551 → 0.2457), and the worst losses are −0.062 and −0.061. The realised leaf count stops varying as well — it moved in 18 of the 27 cells and is now constant, so two runs summarise at the same resolution. Nor is it free: on the engine itself, A-B-A-B on one build, the fit runs 0.83-1.33x across four shapes (200k x 20 1.19x, 200k x 128 0.83x, 100k x 784 1.06x, 20k x 784 1.33x) — a wash to modestly slower. A Python prototype read 0.66-0.96x and was misleading, because it materialised the reordered rows and so bought back sequential access at the price of a second copy of the input. benches/canonical_order.rs splits the cost three ways and shows the mechanism is two opposing effects, not the "coherent inserts split less" first published here: a coherent stream re-descends the same subtree (cache-friendly, down to 0.58x) but fills the leaf budget in one region before the next arrives, so it rebuilds more often -- 3 -> 30 at 50k x 784, max_leaves=8000, which is what makes that shape 1.56x. The key is 3-20 % of the fit and now runs through the SIMD dot kernel; unvectorised it was about twice that. Low-discrepancy walks over the sorted order (van der Corput, round-robin stride) were measured and rejected: 5/16 and 7/16 non-negative against this scheme's 12/16, for an extra magic constant.

Default off, because it relabels. One scoping rule: it applies to a fit / fit_predict that sees the whole dataset — a partial_fit stream never does, so it ignores the flag rather than raising, the same asymmetry as leaf_refit.

Sparse CSR input honours it on the same terms. The key is computed over the CSR rows without densifying, treating a missing column as an explicit zero, so a scipy.sparse matrix and its dense equivalent are sorted into the same order and give the same labels. The merge walk that breaks a key collision needs sorted, unique column indices — scipy's canonical form, but not an invariant of the format — so a matrix that is not in canonical form now raises ValueError naming X.sort_indices(). Silently ordering by a key computed from mis-read rows was the alternative, and it would have looked like a working guarantee. The sparse= leader path (fit_predict_sparse) is a different algorithm and takes no canonical_order.

• leaf_refit= — Lloyd at the micro-cluster level, and an honest account of what it does and does not fix. A leaf CF in every BIRCH-family tree is an absorption history: the rows it swallowed as they arrived, weighted by nothing but arrival order. The tree those rows finished in routes a different partition, and at real compression the two disagree badly enough that row order moves the answer more than the seed does. leaf_refit=k runs k passes of: route every row through the finished tree, rebuild each leaf CF from exactly the rows it wins, drop the entries that win none. Available on Betula, on the free fit_predict, and as CFTree::refit_leaves in Rust. It needs the rows, so a partial_fit stream — which keeps a tree and not the data — ignores it rather than raising.

It is a quality lever at heavy compression, not a fix for order dependence, and the measurements say so. ARI against the labels, medians of seeds 0/1/2, on raw features with everything else at its default: digits at max_leaves=90, kmeans 0.583 → 0.607, ward 0.533 → 0.530, gmm 0.594 → 0.610; MNIST at max_leaves=200, kmeans 0.233 → 0.290, ward 0.216 → 0.258, gmm 0.156 → 0.211. (Re-measured 2026-09-05, replacing numbers that had no single stated setup and one dataset that had been read through the column-major bug above; docs/USAGE.md carries the full table, including what the pass is worth on top of canonical_order.) But the pairwise ARI between two row orders stays at 0.4–0.55 while two seeds stay near 1.0, with or without the pass — the residual order dependence lives in where the prototypes end up, not in what each one summarises, and moving a prototype needs batched routing against a frozen snapshot rather than a better CF for its current position.

Default 0, off, because it relabels: this is a Type-1 door. Two properties are worth knowing before turning it on. Passes are not monotone past 2–3 — each pass re-routes against the tree the previous pass produced, and the worst-leaf gap on the 4000 × 8 fixture runs 2.25 → 1.55 → 0.89 → 0.44 → 1.63 over k = 0..4. And a leaf budget above n is not a free no-op: the partition is the one the tree's greedy descent produces, descent is not the inverse of absorption, and 600 rows at one row per leaf come back as 436 leaves carrying the same total weight.

• KPrototypes grows a third block: directional=. Huang's k-prototypes carries a numeric (n, μ, S) and a category-count histogram per attribute. Some columns are neither — a cyclic hour written as (cos, sin), a heading, a normalised embedding — and this release gives them their own block, summarised by the resultant R = Σ w_i u_i over rows the boundary L2-normalises once each. The dissimilarity becomes ‖Δnum‖² + γ_cat · Σ[x_j ≠ mode_j] + γ_dir · (2 − 2 uᵀc), with c = R/‖R‖. Unit vectors have no mean — two rows a half-turn apart average to zero rather than to a point between them — so the resultant, not a centroid, is the summary; and because ‖u‖ = 1, the block's cost over a whole micro-cluster is 2n_i − 2⟨R_i, c⟩, affine in the summary with no scatter term, the same structure the Lorentzian centroid has in clustering::hyperbolic and what lets a leaf carry the block in one vector. New accessor cluster_directions_ (k × n_directional unit rows). Rust: clustering::{MixedSchema, BlockWeights, MixedRows} join MixedCf; summarize_mixed, nearest_micro and kprototypes now take the schema and a BlockWeights pair instead of a bare γ. Validation moved with it: a numeric column is still always required, but the second block may now be categorical or directional; a column may not be claimed by two blocks; and a one-column directional block is rejected outright, because a direction in one dimension is a sign and a sign is a category. A row whose directional part is all zeros keeps the zero vector — its term is then equal for every prototype, so it stops voting instead of voting arbitrarily, and the same holds for a cluster whose members' directions cancel.
• What the third block buys, measured: the per-row normalisation, and nothing else. The honest baseline is not "no direction" — it is (cos θ, sin θ) written into two numeric columns, and on unit rows ‖u − c‖² = 2 − 2uᵀc exactly, with column scaling by s playing the role of γ_dir = s². Given that same knob on the same grid, the two tie to three digits (0.999 vs 0.999) on unit rows. On rows carrying an amplitude the block reads 0.964 against 0.693 for the naive numeric encoding — and 0.964 again for the user who normalises those two columns by hand. Binning the angle into 8 categories reads 0.509. So directional= is that informed user's preprocessing moved to the boundary, plus a named weight and a unit prototype you can read as an angle; it is not a better model, and on unit rows it provably cannot be. The second argument for the block — that a numeric prototype's +‖μ‖² makes diffuse clusters cheaper and lets them steal members — is real and did not pay: on two adjacent groups at unequal concentration the numeric encoding is 0.011 ahead at both κ = 8 and κ = 2, because for two von Mises–Fisher components of unequal concentration the Bayes boundary is not the angular midpoint and +‖μ‖² moves towards it. Tables and the fixture in bench/RESULTS.md.
• What actually breaks the naive numeric encoding is the amplitude's tail, not its spread. examples/07_mixed_data_kprototypes now carries the third block, on a login-time vector whose angle is the hour of day and whose length is a lognormal session count independent of the segment. Sweeping that length's σ over 0.6 / 0.9 / 1.2 (ARI against the true segments, median of seeds 7/8/9): the raw-numeric encoding reads 0.997 / 0.995 / 0.007 while directional=[4, 5] holds 0.995 flat, because normalisation removes σ before the distance sees it. The collapse is not a spread threshold — at σ = 0.9 the login radius already varies more than monthly_spend does (48 against 29) and costs nothing. It is the tail: at σ = 1.2 squared error is dominated by the longest vectors and the fit spends two of its three clusters on 115 outliers, leaving 1085 of 1200 rows in the third. Normalising the rows by hand reproduces the drop-the-columns fit to three decimals (0.953 at every σ), because unit columns then weigh 1 against numeric columns whose spread is tens of units — which is the scale gamma_dir names and defaults.
• gamma_dir has no literature, and the default says so. γ_cat keeps Huang's ½·mean σ. γ_dir defaults to the mean numeric variance, so that one unit of ‖u − c‖² — which runs over [0, 4] whatever the data — costs one numeric variance. That is a scale-matching convention, not a result; both weights are exposed and both should be swept.

• The k-prototypes++ seeding weight had one mutant nothing could see, and now it has a test. Re-measuring src/clustering/kprototypes.rs whole after the rewrite (162 mutants, 143 caught, 15 unviable, 0 timeout) turned up w[i] * d2[i] → w[i] + d2[i] as a new survivor. The product is what guarantees that a micro sitting exactly on an already-chosen prototype is never chosen again, whatever its mass; the sum drops that guarantee, and the existing spread-out fixture could not tell them apart because there every weight is 1 against a D² of 10⁴. The new fixture puts the mass on the coincident site instead and fails on its first seed under the mutation. The other three entries in mutants-baseline.txt carry over at their new line numbers with their arguments intact.

• method="hyperbolic" — k-means on the Lorentz model of hyperbolic space, for embeddings of hierarchies (Law, Liao, Snavely & Dhillon, ICML 2019). Rows are (d+1) Lorentz coordinates of H^d = {⟨x,x⟩_L = −1, x₀ > 0}; the boundary recomputes x₀ = √(1 + ‖s‖²) per row before insertion, as the directional heads L2-normalize theirs. The head clusters under the squared Lorentzian distance d_L²(x,y) = −2 − 2⟨x,y⟩_L rather than the geodesic arccosh(−⟨x,y⟩_L), because that is the one with a closed-form Fréchet mean — the normalised sum μ = R/|R|_L, R = Σ n_i x_i — and a sum is what a cluster feature carries. d_L² is affine in each argument, so a leaf enters the objective only through (n_i, R_i) and there is no scatter term at all: feature="spherical" reaches the identical partition and identical cost as "full" (measured, ARI 0.6911 to four digits for spherical / diagonal / full), which is true of no other head in this crate. Assignment is argmax_c ⟨R_i, c⟩_L, the cost is Σ_c 2(|R_c|_L − W_c), and the merge increase ΔS = 2(|R_a+R_b|_L − |R_a|_L − |R_b|_L) ≥ 0 is the hyperbolic analogue of Ward's. All of it verified in Maxima before any of it was written (local/scratch/lorentz_identities.mac): the lift lands on the sheet, Lagrange stationarity holds componentwise, S = 2(|R|_L − W) has residual exactly zero, and S ≥ 0, ΔS ≥ 0 survive 400 random weighted sums. Rust: clustering::{hyperbolic_kmeans, HyperbolicKMeans, merge_increase, project_to_sheet, f64_working_radius}. n_clusters=0 is not automatic — the cost falls monotonically in k, as total deviation does for kmedoids. No predict_proba (no density, and the Euclidean softmax fallback would disagree with the head's own labels), no projection= (a code vector has no time-like coordinate), no sparse input (the sheet projection fills column 0 of every row), no refine (a Lloyd sweep is the Euclidean update).

• What the hyperbolic head buys is invariance, not ARI — and the measurement says so. Fixture: a 4-ary depth-4 tree laid out in geodesic polar coordinates, 15 360 points, 16 true groups, threshold=0, max_leaves=2000, median of seeds 0/1/2. On the Lorentz coordinates it reads 0.731 against kmeans 0.407, gmm 0.223, gmm-full 0.663, ward 0.389. But converting the same points to the Poincaré ball and running gmm-full reads 0.817 — the ball is a bounded chart that keeps angle and compresses radius by tanh(r/2), which is what a Euclidean head wants, and on a centred embedding it wins outright. What separates them is a Lorentz boost, an isometry of H^d: the same points, a different origin. At rapidity 3 the ball route falls 0.817 → 0.311 and the tangent-space route 0.767 → 0.375, while this head holds 0.731 → 0.596. The residual 0.14 is the tree's, not the head's — routing and absorption are Euclidean — and raising the budget to one leaf per point makes the head exactly boost-invariant, ARI 0.772 at every rapidity tried. tests/equivariance.rs asserts the exact invariance over a fixed leaf set, with an ambient shift as the control. The f64 working radius is ½ln(2/ε) ≈ 18.4 and is published with the measurement behind it (local/scratch/lorentz_precision.py): relative error of d_L² is 1e−8 at r = 10, 2e−3 at 15, and the sign flips at 18. Tables in docs/USAGE.md.

• method="mfa" — a mixture of factor analysers, mppca with the isotropic residual relaxed to a per-dimension one (Ghahramani & Hinton 1996). Σ_c = W_c W_cᵀ + diag(ψ_c) at rank = q: d − 1 more parameters per component than mppca, against d²/2 for gmm-full. The M-step is Ghahramani & Hinton's, rewritten so the weighted scatter S_c enters only through the d×q product S_c U_c and its diagonal — both of which a leaf answers directly, the second through the new SecondMoment::add_diagonal_scaled, so no d×d matrix is ever formed. Their published S βᵀ (I − βW + βSβᵀ)⁻¹ collapses to W_new = (S U)(M + G)⁻¹ M with G = Uᵀ S U and M + G ⪰ I, which factors by Cholesky where the assembled form needs a general inverse. Nine identities, with the ψ_j left as free symbols, were verified in Maxima at two shapes before any of it was written (local/scratch/mfa_identities.mac, mfa_mstep_collapse.mac; residual exactly zero). rank=0 is a diagonal Gaussian mixture bit for bit. Like gmm and unlike mppca the head is not rotation-equivariant — diag(ψ) is a statement about the coordinate axes — which tests/equivariance.rs now asserts as a property rather than omitting.

• The mfa / mppca envelope, measured in both directions. The two heads dissociate on fixtures built to separate them: mfa reads a separation on a quiet axis that two loud nuisance axes drown (ARI 1.00 against mppca's 0.04–0.34), and mppca reads three lines that differ only in orientation (1.00 against mfa's ≈ 0.00, because ψ absorbs an elongation that belonged in W). On real tables, median of seeds 0/1/2 at rank=2, mfa is behind mppca at full leaf resolution everywhere it was tried — digits 0.562 vs 0.738, MNIST-20k 0.277 vs 0.365, covtype-20k in raw units 0.062 vs 0.087 — at indistinguishable cost (30.9 s vs 37.6 s on MNIST). What it buys is a floor rather than a ceiling: on standardised covtype it converges onto the diagonal head it contains (0.077 / 0.089, gmm's partition exactly) where mppca scores less than half that (0.030 / 0.045). Published in docs/USAGE.md, including why the identical covtype labels are not evidence of a collapsed W: at 1899 leaves in 54 dimensions every one of the 20 000 posteriors is saturated to 1 within 1e-12 for mfa, gmm and gmm-full alike, so agreement there restates the partition, not the covariance.

• method="watson" — a mixture of Watson distributions, for directional data whose sign is arbitrary (Watson 1965; Mardia & Jupp 2000 §9.4.1). p(x | μ, κ) ∝ exp(κ (μᵀx)²) is invariant under x ↦ −x, so eigenvectors, SVD/PCA axes, line orientations and any feature a pipeline signs arbitrarily are clustered as axes, not directions. movmf cannot: handed two equally populated poles of one axis it either spends two components on them or watches its resultant Σ n_i μ_i cancel to nothing. The sufficient statistic is the second moment about the origin, which a cluster feature carries exactly — E[x xᵀ] = Σ_i + μ_i μ_iᵀ — so the E-step term E_{x ∈ leaf}[(μ_cᵀx)²] = μ_cᵀ (Σ_i + μ_i μ_iᵀ) μ_c is closed form and the within-leaf spread is integrated rather than dropped, the same status the Gaussian heads' E-step has. κ < 0 fits girdle (equatorial) components alongside bipolar ones, chosen by likelihood rather than by a flag. Input is auto-L2-normalized as for vmf / spherical-kmeans; n_clusters=0 selects k by BIC. Rust: clustering::{watson, watson_auto, Watson}. Dense path only — a sparse leaf summary is spherical and carries no off-diagonal scatter, so fit_predict_sparse does not accept it.

• What the axial head buys, and what it costs. N = 6 000, max_leaves=300, feature="full", median of seeds 0/1/2, ARI — on an axial fixture with half of each cluster's points at each pole of its axis: 8-D 0.870, 32-D 0.953, 64-D 0.961, against vmf 0.208 / 0.218 / 0.221, gmm-full 0.309 / 0.316 / 0.200 and kmeans 0.233 / 0.244 / 0.232. Its BIC answers k = 4 where vmf answers 8 — one component per pole, the failure the head exists to fix. On digits-PCA20 with a random half of the rows multiplied by −1 and the labels untouched it reads 0.498 → 0.503, while vmf falls 0.631 → 0.173, gmm-full 0.685 → 0.284 and kmeans 0.665 → 0.221. The trade is explicit and ships in the docs: where the sign does carry information the head is a loss (one pole per axis, 32-D: 0.952 against vmf's 0.976), on un-flipped digits it is the worst of the five, a mixture of girdles is not recovered (0.100 on a three-great-circle fixture, where the circles intersect and no axial model can separate them), and it is the most expensive directional head at 5–8× vmf. feature="full", but the margin is narrower than the derivation suggests — measured across max_leaves 4…300, spherical / diagonal / full tie at ARI ≈ 0.95 once every leaf mean has unit norm (~20 leaves up), and full leads by 0.10–0.15 only below that, where a leaf straddles both poles and its mean cancels; fd is the one to avoid at a coarse budget. Tables in docs/USAGE.md.

• The two special-function identities the head rests on were verified in Maxima before it was written. M'(a,b,z) = (a/b) M(a+1,b+1,z) gives the concentration equation g(κ) = (1/d) M(3/2, d/2+1, κ) / M(1/2, d/2, κ) = r̄, and Kummer's transformation M(a,b,z) = e^z M(b−a,b,−z) is what makes κ < 0 computable at all — the ascending series alternates for a negative argument and loses every digit to cancellation, while for a positive one every term is positive and the sum is exact to roundoff. Both residuals are exactly 0 as series identities in z for symbolic a, b. The leading asymptotic M(a,b,κ) ≈ Γ(b)/Γ(a) e^κ κ^(a−b) is not used: against the Maxima table it is 42 % low at d = 50, κ = 50 and still 1.2 % low at κ = 1000, because the expansion parameter at fixed (a,b) is (b−a)²/z ≈ d²/4z, which is not small at moderate d. The series therefore runs everywhere and κ is capped at 10⁴ instead, where its length is bounded. g is strictly increasing with g(0) = 1/d, g(−∞) = 0, g(+∞) = 1, so r̄ ≷ 1/d fixes the sign of κ before any solving and the solver cannot pick the wrong branch. Unit tests check log M and g against 24 Maxima goldens spanning d ∈ [2, 200] and both signs of κ.

• method="dc-median" / method="dc-center" — exact k-median and k-center in the density-connectivity ultrametric (Beer, Draganov, Hohma, Jahn, Frey & Assent, KDD 2023). The same mutual-reachability spanning tree method="hdbscan" takes its hierarchy from, cut for a k the caller names instead of a min_cluster_size. dc(a, b) is the heaviest edge on the path between a and b, so it is an ultrametric and both objectives are solvable, not approximable: k-center by deleting the k − 1 heaviest edges, and k-median by an O(m·k) knapsack over the dendrogram, which exists because a leaf's service cost depends only on which subtree holds its nearest centre and never on which centre it is. Draganov et al. (NeurIPS 2025) run the same recursion for every k at once; this ships the single-k form. Both are asserted against brute force over every C(m, k) centre set below m = 12, scored on an independently computed dc-dist matrix. Neither head emits -1 — they partition, and hdbscan remains the head that answers noise over the same tree. Neither takes n_clusters=0: both costs fall monotonically in k to 0 at k = m, the same reason kmedoids has no automatic mode. Rust: clustering::{dc_clustering, DcObjective}.

• dc-median is a win on the density fixtures; dc-center is mass-blind and ships saying so. ARI / seconds at N = 6 000, max_leaves=600, min_samples=5, median of seeds 0/1/2, k given: moons-0.06 1.000 / 0.021 (spectral 1.000 / 0.203, ward 0.434, kmeans 0.255); moons-0.10 0.889 / 0.014 (spectral 0.691 / 0.383, ward 0.420, hdbscan 0.015); circles 1.000 / 0.015; digits-PCA20 0.725 / 0.020 (ward 0.721, spectral 0.703, kmeans 0.669, hdbscan 0.458); blobs + 5 % uniform noise 0.725 against kmeans's 0.878 — an ultrametric has no notion of a point lying between clusters. dc-center scores 0.000 on moons-0.10, on the noise fixture and on digits: a maximum cannot see a weight, and on a summary an outlier is a low-mass leaf, so at k = 6 its clusters hold [6297, 2, 1, 0] rows against dc-median's [1670, 1240, 1039, 1022]. This is the mass-based answer the CURE probe went looking for and could not find in shrinkage — a different objective, not a repair of the geometric one. The same property read from the other side makes dc-center insensitive to over-specifying k (1.000 / 1.000 / 0.999 / 0.999 / 0.997 at k = 2/3/4/6/10 on moons-0.06, against dc-median's 1.000 / 0.754 / 0.524 / 0.465 / 0.285), because its extra clusters are singletons. Cost is the spanning tree's and within 2 % of hdbscan's own: 0.014 / 0.047 / 0.350 / 1.528 s at max_leaves 300 / 1000 / 3000 / 6000.

• reachability() — an OPTICS reachability plot over the leaves, as the density diagnostic (Ankerst et al. 1999). Returns a ReachabilityPlot with order / reachability / core_distances / weights and a labels_at(eps) that reads the DBSCAN* cut off the plot. It is the method="hdbscan" hierarchy, not a second opinion about it: OPTICS with no ε cutoff is Prim's algorithm on the mutual-reachability graph, so the sweep walks the same spanning tree that head cuts, and the crate builds it once — clustering::hdbscan::mutual_reachability is now shared by both. Two theorems carry it, asserted rather than measured: reachability[1..] is a permutation of that tree's edge weights, and cutting the plot at ε gives exactly its components under ε. That is also why the reachability is the mutual max(core(p), core(q), d(p,q)) rather than Ankerst's asymmetric max(core(q), d(q,p)) — the asymmetric form pictures a different tree. One sweep position per leaf, so a valley's width is a leaf count; weights is published next to it for exactly that reason. Rust: clustering::{optics, Reachability}.

• The diagnostic's cost is the leaf count, and it loses the fixture that is about noise. Median of seeds 0/1/2. Cost does not see N — 0.0024 / 0.0027 / 0.0027 / 0.0028 s at 5 k / 20 k / 80 k / 320 k rows against a 0.005 / 0.007 / 0.017 / 0.054 s fit — but it does see the leaf budget quadratically: 0.0003 / 0.0028 / 0.0284 / 0.2970 s at 99 / 290 / 926 / 2883 leaves. graph_degree bounds that and is not free: at 2883 leaves graph_degree=16 runs 8× faster but its cut agrees with the exact plot's at only ARI 0.518, and 32 / 64 / 128 all sit at 0.923 — an approximate graph is a plot of a slightly different tree, not a cheaper plot of the same one. Read as a partition against sklearn.cluster.OPTICS on the raw points (N = 6 000, max_leaves=600, best ε on each side's own grid): blobs 0.452 vs 0.448, moons 0.997 vs 0.978, circles 0.997 vs 0.978, blobs + 5 % uniform noise 0.687 vs 0.753, at 0.024–0.041 s against 4.9–6.3 s. The compression is not what costs the noise fixture — raising the budget to one leaf per point takes it only to 0.707, so the residual is the convention: DBSCAN* has no border points, and the mutual reachability is a stricter link. What the budget does fix is the mass-weighted core distance — at 274 leaves 41 % of the true noise points sit in a leaf with core distance 0, which can never be labelled noise; at 3899 leaves that is 2 %, and the ARI moves 0.518 → 0.686.

• method="fuzzy-cmeans" — weighted fuzzy c-means over the leaf centroids (Bezdek 1981). The only soft head in the crate that fits no density: predict_proba returns the memberships u_ij ∝ d_ij^(−1/(m−1)), a partition of unity over the centres rather than a posterior. Exact on the summary by the same identity kmedoids uses — tying the membership within a leaf turns Σ_{x ∈ leaf i} u_ij^m‖x − c_j‖² into u_ij^m · n_i · d_ij with d_ij = ‖μ_i − c_j‖² + S_i/n_i, so the leaf-level J_m is the point-level one. The new fuzzifier keyword is the exponent m > 1 (default 2.0); it is validated rather than clamped, because at m ≤ 1 the exponent 1/(m−1) inverts the rule into "farthest centre wins". n_clusters=0 selects by Xie–Beni, this family's own validity index.

• The fuzzy head's hard partition is a loss, and it ships saying so. ARI against kmeans at max_leaves=4000, feature="spherical", median of seeds 0/1/2: blobs 0.847 vs 0.864, aniso 0.535 vs 0.545, varied 0.531 vs 0.540, highdim 1.000 vs 1.000, digits 0.483 vs 0.467 (all at m = 1.3; m = 2.0 is worse everywhere except varied). The intuition that it should win where clusters overlap is backwards: on 6 000 × 8 blobs at cluster_std 3/4/5 it reads 0.990 / 0.913 / 0.623 against kmeans's 0.990 / 0.922 / 0.804, and the confidence AUC — how well assignment_confidence separates right from wrong labels — falls with it (0.975 / 0.916 / 0.836 vs 0.985 / 0.931 / 0.881), so the membership is not a better boundary detector either. Mechanism: no membership is ever zero, so every point pulls every centre toward the grand mean. Measured on the cluster_std=5.0 fixture, the mean centre-to-centroid distance runs 13.75 for kmeans and 13.19 / 12.98 / 11.88 at m = 1.3 / 2.0 / 3.0 against a true 13.67. Take the head for the membership, not for the labels. Rust: clustering::{fuzzy_cmeans, fuzzy_cmeans_auto}. Tables in docs/USAGE.md.
• method="kmedoids" — weighted k-medoids over the leaf centroids, by eager FasterPAM. The centre of a cluster is one of the summary's own micro-clusters rather than an average: an exemplar you can show, and a centre that stays on the data manifold where a mean need not. Schubert & Rousseeuw 2021; the ΔTD(i, c) = shared(c) + corr_i(c) decomposition prices one candidate against every medoid slot in a single O(m) pass. Exact on the summary under the whole-leaf restriction every shipped head accepts — Σ_{x ∈ leaf i} ‖x − μ_j‖² = S_i + n_i‖μ_i − μ_j‖², so the leaf-level total the swap search minimises is the point-level sum of squares up to the constant Σ_i S_i. Note the square: classical PAM minimises absolute distance, which has no closed form in a cluster feature, so the squared objective is the one a summary can answer exactly.

On digits (feature="spherical", threshold=0.0, ARI median [min–max] over seeds 0–4) it reads 0.554 [0.554–0.570] at one leaf per point against kmeans's 0.467 [0.443–0.571] — better on the median and on the spread, because the centre cannot drift into the gap between two digit classes — and 0.520 [0.468–0.520] against 0.487 [0.381–0.523] at 296 leaves. It loses at a coarse summary, 0.219 against 0.240 at 115 leaves, where there are too few candidate centres to place ten of them well. The O(m²) swap pass costs 6× kmeans at 1797 leaves and nothing measurable at 296. On the synthetic fixtures at max_leaves=4000 the two are level (blobs 0.864/0.864, aniso 0.545/0.548, varied 0.540/0.548, highdim 1.000/1.000, median of seeds 0/1/2).

Two consequences the head states rather than hides. n_clusters=0 is not free auto-k: total deviation falls monotonically in k, so this objective cannot choose one, and the automatic arm runs the medoid silhouette (dyn_msc, Lenssen & Schubert 2024) — a different objective. And refine is a no-op, because a Lloyd sweep is the k-means update and would move each medoid to the mean of its cluster, off the data. fit_predict_sparse rejects the head: on its flat leader summary a row's distance to a single micro-cluster is dominated by that micro-cluster's norm rather than by the overlap that knows the topic, measured at ARI 0.017 (own rule) and 0.002 (micro-cluster route) against 1.000 for kmeans — and 1.000 for kmedoids itself through Betula.fit, which builds a real CF-tree from the same CSR. Rust: clustering::kmedoids → Pam. Tables in docs/USAGE.md.

• DenStream.drift_ and DbStream.drift_ — an ADWIN change detector on the streaming heads. decay is a schedule: it forgets at a fixed rate whether or not anything changed, and a wrong λ is silent in both directions. The detector (Bifet & Gavaldà, SDM 2007) answers the other question — did the stream change — with a stated false-positive ceiling δ = 0.002 rather than a tuned threshold, and reports {alarms, last_alarm, distance, window}. It reports; it does not act: an alarm prunes nothing, promotes nothing and relabels nothing, so no existing labels move. The statistic is the distance from each incoming point to the nearest micro-cluster in units of the micro-cluster radius, which is free (both online steps already compute it) and scale-free (a raw distance stops detecting anything on data scaled down far enough, because ε_cut's additive term is in absolute units). Measured over jump / spread / split / gradual changes on both heads, three seeds: reported in 10 of 10 cells, three quarters of them inside the first 32-point check clock, at a false-positive rate of 0.00050 on stationary streams. The obvious alternative — a novelty indicator, 1 if the point opened a micro-cluster — was implemented first and measured to raise zero alarms in five of six (λ, ε) settings on a 50σ jump: the impulse lasts ~5–10 points and a bounded statistic cannot show a gap that short. New bench/drift.py, tables in bench/RESULTS.md.

• The automatic-k ceiling is no longer one number for every head, and auto_k_max overrides it. n_clusters=0 was bounded at 20 for all of them, so no n_clusters=0 path in the library could return more than 20 clusters. The bound exists because a sweep refits the whole head at every candidate k — Σ_{k≤K} k = O(K²) — but ward, average, weighted, centroid and median build one dendrogram and score its cuts, and xmeans stops on its own split test, so they were paying a cost guard for a cost they do not have. Those six are now bounded only by the leaf count. Measured on 480 leaves in 64 dimensions holding 120 true groups, ward goes from 5.7 ms at k = 2 and ARI 0.009 to 23.8 ms at k = 120 and ARI 1.000. method="ward" / "average" / "weighted" / "centroid" / "median" labels change at n_clusters=0 whenever the data has more than 20 groups.

The sweeps keep the ceiling, because lifting it there costs 8–32× — kmeans 46 ms → 1.4 s, gmm 0.5 s → 4.1 s on the same fixture — and take auto_k_max (0 = the default) as an explicit opt-in. A selection that lands exactly on its ceiling now raises a UserWarning: an argmax on the last candidate is evidence the search stopped early, not evidence about the data. On k-means-shaped data method="xmeans" reaches the same partition as the fully-swept kmeans 120× faster (12.1 ms against 1449 ms) and is the thing to reach for before auto_k_max. Table in docs/USAGE.md.

• method="xmeans" — x-means proper (Pelleg & Moore, ICML 2000), the first auto-k head that can return more than 20 clusters. It alternates k-means with a per-centre split test — each centre is asked separately whether the leaves it owns are better described by two centres than by one, scored by BIC on that subset alone — and stops when no centre wants to split. n_clusters is an upper bound, not a target; 0 bounds it only by the leaf count, because the split test rather than a cost guard is what stops the recursion. The sweep behind method="kmeans" cannot do this: it costs O(k_max²) full k-means over every leaf, which is why AUTO_K_MAX = 20 exists and why no other auto-k head in the library can exceed 20.

The split test has a threshold and it is a function of d. A balanced binary split costs n·ln 2 of mixture-weight likelihood and buys ½·n·d·ln(S₁/S₂), so it is accepted only when the cut captures more than 1 − 2^(−2/d) of the region's sum of squares — half of it at d = 2, 0.13 at d = 10, 0.022 at d = 64. A cut through a cloud that is round at every scale captures about 0.64/d, always less than the 1.39/d required. The recursion therefore starts at k = 2: a greedy splitter has no way back from a refused split, and at k = 1 the whole answer rides on the one comparison that rule answers worst, since a layout of many well-separated groups is itself close to round. From k = 1 the head collapses to a single cluster in five of twenty (k*, seed) cells, in every seed at k* = 60, and on a 3×3 grid of nine equal 2-D blobs; from k = 2 every one of those returns the true count. Pelleg & Moore, ELKI and pyclustering all start at 2. Measured on random blob layouts (4 leaves per blob, median of seeds 0/1/2), |Δk| against the true count with the sweep at its cap of 20: d = 10, k* = 30 → sweep 10.0 (capped) in 4.8 ms, x-means 0.0 in 0.8 ms; d = 64, k* = 30 → sweep 10.0 in 9.0 ms, x-means 0.0 in 1.4 ms. What survives is d = 2 at a large count, where one cut must capture half the scatter: d = 2, k* = 30 → sweep 10.0 (capped), x-means 6.7. Table in docs/USAGE.md.

• Betula.near_duplicate_pairs(..., neighbors=k): a cross-leaf recall pass. The within-leaf scan cannot see a duplicate pair whose two rows were absorbed into different leaves — at 20 000 rows with 400 planted pairs it misses 135 of them outright — and neighbors scores each leaf against its k nearest leaves as well, recovering about 40% of the split pairs at k = 8. 0 (the default) is the previous behaviour unchanged. Table in bench/RESULTS.md.

• Betula accepts max_leaves as a fraction of the row count: a float in (0, 1) resolves to ceil(frac · N) at fit time, alongside the absolute integer cap. (Estimator only — resolving the fraction lives in the Python wrapper, while the free fit_predict / fit_predict_sparse are the engine functions re-exported verbatim and take the integer.) One parameter carries both forms, discriminated by range rather than by a second flag, which is ELKI's -cftree.maxleaves convention (its own default is 0.05) and makes a configuration independent of dataset size. A fraction is undefined for partial_fit, which never sees a final N, and raises there rather than resolving against a batch; memory_budget_mb still overrides either form. Reported by effective_max_leaves_ as before.

• Mapper linkage option link="bhattacharyya": microclusters inside a cover bin are linked on the Bhattacharyya distance between their Gaussian surrogates — the centroid gap divided by the pair's own spread — instead of centroid distance alone. It attacks chaining, where a thin bridge of sparse microclusters merges two dense regions: on a dumbbell fixture at three seeds it crosses the bridge in none of 63 runs while the centroid rule merges the two lobes into one at link_scale = 3. The default is unchanged. The bounded Hellinger distance was measured on the same fixture and refuted — it saturates at the leaf spacing a CF-tree produces, so any link_scale past ≈ 1.1 links the whole bin — and is not exposed. Table in bench/RESULTS.md.

• mixture_w2 and Betula.summary_w2, the mixture-Wasserstein distance MW₂ between two fitted Gaussian mixtures (Delon & Desolneux, SIAM J. Imaging Sci. 13(2), 2020). The pair cost is the closed-form Bures W₂; the transport over the component grid is solved exactly by a transportation simplex, so a mixture scores exactly zero against itself rather than nearly zero. summary_w2 applies it to two models' leaf summaries and is a drift metric that needs no labels, no shared points and not even a common leaf count — on a four-blob drift sweep the label ARI stays at 1.0000 in twelve of thirteen scenarios while MW₂ tracks every one of them, matching √(w·δ² + null²) to within 2.5%. mixture_w2 takes parameters rather than an estimator, so an sklearn or ELKI fit can be compared without conversion. Table in bench/RESULTS.md.

• Diffusion-maps α-normalization on the spectral head — measured and refuted, default unchanged. The head's Laplacian sits over mass-weighted leaves, so sampling density enters it twice, and Coifman & Lafon's A_ij / (q_i^α q_j^α) is the standard correction. Over four fixtures and a density-imbalance sweep up to 1:25, medians of three seeds, α never wins a cell: wherever the columns differ at all, α = 0 ≥ α = 0.5 ≥ α = 1, and on three sheared ribbons ARI falls 0.8892 → 0.7516 → 0.6903. The reason is that the correction is already applied — the affinity is self-tuning (Zelnik-Manor & Perona local scaling), its σ_i is a local density estimate, and α is derived for a fixed-bandwidth kernel, so stacking them divides the same bias out twice and over-corrects. DIFFUSION_ALPHA stays 0.0 and now records that as a measurement.

The affinity's neighbour count was measured on the same fixtures and turns out to be load-bearing: k = 24 collapses two of four, and k = 7 beats the shipped k = 10 on anisotropic shapes (0.9300 against 0.8892) while tying elsewhere. knn_degree now carries the Γ-convergence floor that makes a choice defensible — García Trillos & Slepčev require k_n / log n → ∞, and at the head's SPECTRAL_MAX_NODES = 256 cap the shipped 10 is 1.8 log n while k = 5 (the best cell measured) is below the threshold. The default is unchanged at 10 because knn_affinity is shared with the Leiden head and four synthetic fixtures are not the evidence base for a constant that relabels two heads at once. Tables in bench/RESULTS.md. - gap_statistic: a k-selector that is allowed to answer k = 1. Calinski–Harabasz is undefined at k = 1 and the distance-based indices hallucinate structure in noise — Schubert, Stop using the elbow criterion for k-means (SIGKDD Explorations 25(1), 2023), Table 1 records Silhouette and Davies–Bouldin reporting 3–22 clusters in pure noise where BIC reports one, and the gap statistic doing so only partly. This computes the gap on the leaf summary: pooled within-cluster dispersion from S_i + n_i‖μ_i − c‖² (exact, O(ℓ k d), no second pass over the points) against a uniform null over the data's bounding box that is re-summarized at the same leaf budget, so both sides pay the same quantization error and the gap measures structure rather than compression. The reference draws are shared across k, which pairs the comparison. Measured on that table's two nulls (1000×5, seeds 0/1/2): a single Gaussian answers k = 1 on every seed, uniform noise on two of three — and the third is a near-tie, gap(1) − gap(2) = −0.007 against −2.83 for a real two-cluster fixture, i.e. undecided rather than confidently wrong. Two- and four-blob controls come back at 2 and 4 on every seed. docs/USAGE.md now carries the which-selector-when table (BIC parametric and authoritative, CH cheap but k > 1 only, gap for the k = 1 capability without a distributional assumption). - tree_report() and estimate_threshold: why the tree collapsed, and whether it mattered. Betula.tree_report() answers "why is my tree collapsing?" from the leaf summary alone — budget fill, the heaviest leaf's share of the mass, its width against a typical leaf, and leaf-mass quantiles — and turns them into plain-language findings. The width is the part that is not obvious: mass concentration alone does not distinguish a real loss from a faithful summary, because bench/size_imbalance.py's structured and flat fixtures both put 80% of the mass in one leaf at a 93–96% full budget. The heaviest leaf's radius over the median leaf radius separates them 3×: 0.53 / 0.75 for structured (two clusters merged inside that leaf) against 0.17 / 0.27 for flat (nothing inside to lose), at max_leaves 250 and 4000, medians of seeds 0/1/2. On MNIST the mass trigger fires at exactly the two leaf budgets whose ARI is materially below the best cell (0.185 and 0.235 against 0.284) and stays quiet from ×22 compression down. Passing X adds an A-BIRCH gap-statistic threshold estimate (betula_cluster.estimate_threshold, Lorbeer et al. 2018) next to the threshold in use, so a tree that settled above twice the sampled estimate is named as absorbing points from more than one cluster. The estimate is advisory — max_leaves is still the knob that binds — and it reports each of the paper's assumptions the data breaks (comparable cluster sizes, well-separated clusters) rather than assuming them. It departs from Tibshirani, Walther & Hastie's original selection rule for a measured reason: on four blobs 40σ apart the gap curve is non-monotone (gap(2) = −0.186, gap(3) = −0.274, gap(4) = +3.379), so the "first k that beats k+1 within one standard error" rule stops at 2. Restricting that rule to the k within one standard error of the best gap returns 4, still returns 1 on a single Gaussian blob, and declines to certify uniform noise. - outlier_scores(X, metric="mahalanobis"): whiten by the cluster's covariance, not its trace. The default metric="radius" divides the centroid deviation by one scalar RMS radius, which is the trace of the cluster's pooled covariance — so an elongated cluster's short axis is judged by the length of its long one. The new metric pools the leaves by the parallel-axis theorem (Σ = Σ_l w_l (Σ_l + δ_l δ_lᵀ) / W, ridged at the GMM head's VAR_FLOOR_REL) and whitens against a Cholesky factor of that. The two are calibrated rather than merely both being z-scores: on an isotropic cluster they return the same number (checked on the 2⁵ hypercube corners and in 1-D, to floating point), so the refinement moves a score only where the cluster has a shape. On the DAMI fixtures, medians of seeds 0/1/2: sheared clusters 0.596 → 0.748 ROC-AUC (average precision 0.059 → 0.181), blobs+uniform 0.991 → 0.995, unequal-density unchanged at 1.000 — nothing regresses. The off-diagonal terms are the whole gain: a per-dimension variance, which is what "reuse the second moments the CF carries" first suggests, reaches only 0.642, because a 37°-rotated ribbon's coordinate variances differ by under 2×. Costs O(k d³) once plus O(d²) per row against the scalar path's O(d), so the default is unchanged. find_outliers passes metric through. - outlier_scores measured on the DAMI 2025 axes: where it wins, and the shape it cannot see. Sanchez Vinces, Schubert, Zimek & Cordeiro's Clustering-based outlier detection names the axes — detection quality, resilience to parameter variation, scalability, and auto-filtering by an internal index — and three of the four are now measured against IsolationForest and LOF over three fixtures, medians of seeds 0/1/2. The per-cluster z-score wins the local case decisively: on clusters of unequal density it reaches ROC-AUC 1.000 / average precision 0.998 against IsolationForest's 0.887 / 0.175, with the smallest parameter spread of the three detectors (0.001 over a 12-cell sweep). It loses on anisotropy — ROC-AUC 0.596 on sheared clusters, and identically 0.596 at feature="diagonal" and feature="full", which is the mechanism: the score divides by one scalar RMS radius, so a cluster's short axis is judged by the length of its long one and the covariance the tree already carries is never consulted. A Mahalanobis radius on the same partition lifts that to 0.762, filed as its own task. The parameter cell can also be chosen without labels: taking the largest validity()["calinski_harabasz"] over the sweep lands within 0.033 ROC-AUC of the sweep's best cell on all three fixtures — so the resilience spread is filterable rather than fatal. No behaviour changed; docs/FEATURES.md and bench/RESULTS.md now state both the win and the limit. - fidelity::summary_mmd: a label-free, head-independent fidelity number for the leaf summary. Every other quality number the crate reports needs something the summary cannot supply — ground truth labels, a head, or a value of k — and mean_sq_radius, the one that needs none of them, is the exact k-means fidelity and falls monotonically forever. Maximum mean discrepancy asks the distributional question instead. The leaf surrogate is the Gaussian mixture Σ (n_i/N)·N(μ_i, s_i I) with s_i = S_i/(n_i·d), the same isotropic leaf cloud the scale-space head mollifies with, and the expected Gaussian kernel between two such clouds is (1 + s/h²)^(−d/2)·exp(−‖m‖²/(2(h² + s))), so the MMD against a raw sample needs no sampling at all. Verified against Monte Carlo before use: worst |z| = 3.02 over 24 stochastic cells at 4·10⁶ draws, the 8 deterministic cells exact to 10⁻¹².

Measured, it finds the knee the quantization error hides. On six overlapping Gaussians in 16 dimensions, ARI jumps 0.593 → 0.976 between budgets 16 and 32 and is then flat; summary_mmd falls 0.0665 → 0.0065 across the same steps and is then flat within ±5% for five more budgets, while mean_sq_radius slides smoothly 212 → 80 with no feature at the knee. It is not a monotone axis — in two dimensions it rises from 0.0054 at 8 leaves to 0.0092 at 32 before falling, because coarsening also lets a leaf's ball fit a Gaussian blob better than a two-point leaf's does, and the ARI column dips at exactly the same budget. Exposed as Betula.summary_mmd(X, bandwidth=None), which needs only a built tree — no labels, no k, no head, so a partial_fit tree can be scored. O(m² + m·M + M²) kernel evaluations: a diagnostic to run at a few budgets, not inside a fit. Tables in bench/RESULTS.md. - The symmetry group of every (feature, head) pair is now stated and enforced. A head answers a question about the data, not about the coordinates it arrived in, and which frame changes it may ignore is a property of its model that is easy to lose by accident — to a bandwidth in absolute data units, a covariance floor applied before whitening, or an initialisation that reads coordinate 0. tests/equivariance.rs runs all eleven Euclidean heads in five frames (translation, rotation, uniform scaling, an axis swap, and the identity) and compares the partitions up to relabelling, with the two exceptions asserted to differ rather than skipped: gmm_diagonal is not rotation invariant because an axis-aligned covariance is a claim about the axes, and the directional heads are not translation invariant because a direction is read from the origin. The audit found no head violating its claimed group — reportable only because the suite carries its own controls, including one that superimposes two of the four groups and requires most heads to notice, since blobs separated far enough to be found in any frame would pass an invariance suite without exercising a symmetry.

The measured half is the second axis: the leaf summary is as much of the group as the head is. A symmetry the feature has already discarded cannot be recovered downstream, so gmm_full — a rotation-invariant model — is rotation invariant on Full leaves, not invariant on Diagonal ones (the default cov_dense keeps only diag(variance), and the off-diagonal scatter these clusters live along is gone), and invariant again on Spherical ones for the trivial reason that an isotropic summary has nothing a rotation can change. Choosing feature="diagonal" for memory is therefore also a choice about which frames the answer may depend on.

Auditing gmm_toeplitz — neither Euclidean nor directional, since an AR(w) covariance says the coordinates are an evenly spaced sequence — needs a third kind of frame change: shifting or rescaling the whole sequence leaves that model alone, while reordering the axes is a different model. It is exhibited on two AR processes differing in nothing but their autocovariance, each window standardised so the marginal mean and variance carry no signal, which the head separates exactly; sending lag-1 neighbours 64 apart destroys the band structure and the answer with it. Tables in docs/MATH.md. Tests only — no behaviour changes. - dyn_msc: medoid-silhouette clustering that chooses its own k. (Rust API — clustering::dyn_msc, clustering::MedoidClustering.) Lenssen & Schubert's DynMSC (Inf. Syst. 120, 2024): optimise the medoid silhouette by swaps, then sweep k downward reusing the medoid set and keep the best-scoring level. Its objective is exactly the already-published validity::medoid_silhouette, asserted rather than assumed, so the head cannot win on a number the metric would score differently.

It is the crate's second automatic k, and it does not replace the first. Over 24 measured cells it is 11 wins to Calinski–Harabasz's 9 with 4 ties, and the split is by cluster shape: on round clusters DynMSC recovers k in 9 of 12 cells and gives the better partition whenever both agree on k, while on clusters stretched 4× along one axis CH returns the top of the k range it was handed in 5 of 12 cells — the monotone failure of a variance ratio on non-spherical data — where DynMSC's error is bounded by collapsing toward 2. Cost scales as m^2.4–m^2.6, about 0.37 s at 960 leaves, so it belongs to the same hundreds-of-leaves envelope as the other quadratic heads. ward_hac_auto keeps the default and there is no Python surface. Tables in bench/RESULTS.md.

• Selection: PLSCAN's persistence-chosen minimum cluster size, alongside excess of mass. (Rust API — clustering::hdbscan_selected, clustering::Selection.) Raising HDBSCAN*'s min_cluster_size never moves a merge, it only prunes branches that fail to reach the size, so one dendrogram already contains every clustering the parameter can produce. Each segment therefore has a size interval over which it is a leaf cluster; summing those lifetimes rates a candidate size, and Selection::Persistence reports the clustering at the argmax (Bot, McInnes & Aerts, arXiv:2512.16558). Hdbscan gained selected_size so the chosen number is visible.

Measured over 50 cells, it is a conditional win and ships with the condition stated: on clusters that differ in density it wins 14 / loses 0 / ties 11 and its ARI spread across min_samples is never larger; on clusters that differ in size it wins 3 and loses 22. The excess-of-mass default is unchanged and there is no Python surface — a method= string cannot carry a condition. Part of the paper's benefit is already bought by the CF tree, since a leaf of w points cannot express a cluster finer than w; the two arms converge as leaves get heavier. Tables in bench/RESULTS.md.

• dendrogram_below: exact agglomerative clustering below a height bound, with a certificate. (Rust API — clustering::dendrogram_below, certificate_radius, BoundedDendrogram.) A cut at k clusters only needs the merges below one height, and those can be computed over a candidate radius graph instead of all pairs — exactly, not approximately. The argument is one inequality: the minimum over cross pairs is at most their mean, ‖Δμ‖² + S_A/n_A + S_B/n_B. Average and weighted linkage get r² = h because their value is that mean; Ward gets r² = h/w_min because its value carries a mass factor and a cluster's scatter is the sum of the D4 of the merges inside it. Centroid and median get no radius at all — their value is ‖Δμ‖², which bounds neither spread, so two concentric shells merge low with no close cross pair anywhere. That case returns CertificateError::NoRadius rather than a plausible wrong answer, and a cut above the certified height returns None rather than one served from a truncated dendrogram.

Measured at 7.9× the dense driver on 32 000 leaves (10.3 s against 81.9 s, median of 3), 7.6× at 8000 and 5.9× at 2000. The graph is still an O(m²·d) scan — 26–42 % of the total — because an exact radius query has no sublinear structure here, and an approximate k-NN index would forfeit the certificate for the sake of the one edge it drops. Full tables, including the numeric cross-check that the three bounds are attained, in bench/RESULTS.md.

• BregmanBetula: the CF-tree in a Bregman geometry, as a second estimator. Lane A shipped the Rust side in this release — BregmanCf, D4_φ, Bregman k-means / Ward / EM — with no Python surface, because the wiring is a type-1 door. It is now decided and recorded in docs/adr/004-bregman-public-api.md: a separate class, not a divergence= keyword on Betula.

The reason is that feature= names a covariance model and every one of its values summarises the same Euclidean geometry, while a Bregman feature names the geometry. Collapsing the two axes makes meaningless combinations writable — feature="kl" with method="gmm" hands a Gaussian head a Bregman information to read as a variance, absorb="chi2" applies a variance prior to a quantity that is not one. With a second class they are not rejected; they cannot be typed.

BregmanBetula(divergence=…, method=…, beta=…) takes divergence in {"euclidean", "kl", "itakura-saito", "logistic"} and method in {"kmeans", "ward", "mixture"}. beta is the mixture's inverse dispersion in nats of divergence, and it is rejected rather than ignored under the other two heads. The domain each φ needs — x > 0 for KL and Itakura–Saito, x ∈ (0, 1) for logistic — is validated at the Python boundary and the error names the offending row and column, because BregmanCf::push only debug_assert!s it and a release build would otherwise return NaN. divergence="euclidean" is kept deliberately: squared Euclidean is a Bregman divergence, so it reduces to the shipped path, which is what makes the agreement between the two estimators a test rather than a claim.

• balance: a mass-balanced leaf budget, for the data where the geometric one fails. The absorption gate is purely geometric — one global radius, raised by the rebuild until the leaf count fits — and a single radius cannot serve two densities. Once it passes a dense region's diameter that region collapses into one leaf in a single step while sparse regions keep splitting: measured at 80 000 of 100 000 points in one leaf at every budget from 250 to 4000, with the budget itself 90–96 % filled. The budget was never under-used; it was spent on whichever points happened to be far apart.

balance = b caps a leaf at b × (total mass / max_leaves). The cap is enforced at both sites that combine mass, since either alone undoes the other — absorption refuses a full entry and starts a new one, compaction skips a pair that would overflow — and max_leaves stays a hard bound: a rebuild that cannot reach its target under the cap merges over it rather than leave the tree over budget. The ideal is read from the root CF, so a stream tightens it as data arrives instead of needing n up front.

On the fixture that motivated it — two true clusters inside a dense core, five diffuse minorities around it — kmeans goes from ARI 0.4174 to 1.0000 at every budget from 250 to 4000, matching sklearn's k-means on the raw points, while the realised leaf count barely moves and the control fixture does not regress.

Default off, because it is a lever and not a free win. Across digits, covtype-20k and mnist-10k × three heads × three budgets, 19 of 27 cells improve — but the useful statistic is the heaviest leaf's share of the mass, not the average: every cell with top1 ≥ 0.5 gains +0.08 to +0.25 with all three heads agreeing, and every cell with top1 < 0.1 moves within seed noise in both directions. max(microcluster_weights_) / sum(...) is the one-line diagnostic that says which case you are in. Full table in bench/RESULTS.md. - bregman: the cluster feature, generalised from squared Euclidean to any Bregman divergence. For d_φ(x,y) = φ(x) − φ(y) − ⟨∇φ(y), x−y⟩ the within-cluster Bregman information S_φ = Σ w_i d_φ(x_i, μ) obeys exactly the identities the Euclidean scatter does — Maxima residual 0 for a generic differentiable φ — and the mean merge is unchanged, which is Banerjee's characterisation rather than a coincidence: the arithmetic mean is the right-sided Bregman centroid for every φ. So (n, μ, S_φ) is a commutative monoid for every φ and the tree machinery is untouched.

The trait's shape is dictated by a measurement. Σ w φ(x) − W φ(μ) is the exact analogue of BIRCH's SS − n μ² and cancels the same way — but the recurrence alone does not save you either, since both of its new terms evaluate d_φ at arguments that nearly coincide, and expanding d_φ from its definition subtracts two nearly equal values of φ one level down. On the logistic divergence at a tight cluster the naive form loses everything (relative error 1) and the expanded recurrence diverges, returning an information 10¹¹ times too large. So BregmanDivergence requires a hand-written cancellation-free divergence() whose contract is accuracy when x ≈ y, and keeps phi() as a test oracle for well-separated arguments only.

Ships SquaredEuclidean, KullbackLeibler, ItakuraSaito and Logistic; BregmanCf implementing ClusterFeature; and the two measures a tree needs — BregmanCentroid (D0_φ, routing) and BregmanIncrease (D4_φ, absorption). Both live in bregman rather than distance because, unlike everything there, they are not feature-agnostic: they have to know which φ built the feature. The regression test that matters is that φ(t) = t² reproduces the shipped Spherical feature and that D4_φ collapses to VarianceIncrease — the generalisation must not have changed the Euclidean answer.

Two Phase-3 heads over that summary: clustering::bregman_kmeans and clustering::bregman_agglomerative. Both are exact on the leaf summary for the same reason the Euclidean ones are -- the bias-variance identity holds for every φ, so a leaf is scored against a candidate centre without revisiting a point -- and the k-means seeding uses the whole leaf potential S_i + n_i·d_φ(μ_i, c) rather than the centroid's share alone. Two things deliberately do not carry over: d_φ is not a metric, so the Hamerly bounds that accelerate the Euclidean k-means have no analogue and the Lloyd loop is the plain one; and Bregman-Ward is not reducible, so the agglomerative head runs Anderberg (see ADR 002 below). At φ(t) = t² the agglomerative head reproduces the shipped Euclidean Ward labels exactly, at every k tested.

The Anderberg driver added in 0.7.0 for the non-reducible Euclidean linkages is now node-agnostic and shared by both, rather than duplicated -- a behaviour-preserving change, pinned by the existing SciPy and Lance-Williams cross-checks.

clustering::bregman_em closes the family: soft Bregman clustering is EM for the matching regular exponential family (Banerjee's bijection), so it is the diagonal-GMM skeleton with d_φ where the Mahalanobis term was. Two things about it are worth stating.

E_{x∈leaf}[d_φ(x, μ_k)] = S_i/n_i + d_φ(μ_i, μ_k) is exact for every φ, so the expected complete-data log-likelihood needs no approximation of the within-leaf shape at all; tying the responsibilities within a leaf is the only approximation, which makes this exact variational EM with a monotone lower bound. A consequence that is testable and surprising: S_i/n_i is the same for every k, so given the centres a leaf's internal spread cannot move its responsibilities — only the value of the bound. (The fitted result still depends on S, through the k-means++ warm start.)

And beta, the shared inverse dispersion, is a parameter rather than a fitted quantity — fitting it needs the family's log-partition function, which has no form generic in φ. It is also what decides whether the head can separate anything, because at fixed dispersion separation is measured in nats of divergence, not in coordinates. Itakura–Saito is scale-invariant, so three groups whose centres look far apart can be 0.33 nats apart, and a β = 1 mixture correctly reports that they overlap and collapses all three means onto the global mean. That is the model being honest rather than the optimiser failing, and there is a test pinning both halves of it.

Rust-only for now; nothing is wired into Method or the Python API yet. - ADR 002 — a Bregman-Ward HAC must use Anderberg, not the nearest-neighbour chain. NN-chain is exact only for a reducible linkage, and generalising Ward's criterion to a Bregman divergence does not preserve reducibility. Measured rather than assumed: zero violations at d = 1 across 213 million admissible triples and 16 different φ, and violations from d = 2 for every φ except squared Euclidean. Inside every pooled counterexample, no coordinate ever violated its own 1-D reducibility — the failure is purely that reducibility is not preserved under addition. φ(t) = tᵖ is clean at p = 1.99 and p = 2.0 and already violates at p = 2.01, so there is no safe neighbourhood of t².

It costs a wrong dendrogram, not a wrong triple: against exact Anderberg over 3 000 instances per cell, squared Euclidean differs in 0 of 9 cells, while Itakura–Saito at d = 20, m = 12 differs in 1.0 % and the exponential divergence in 1.2 % — a rate that grows with m, on a head that runs at m in the thousands — and when it fires the answer is destroyed rather than perturbed (ARI 0.10 at k = 4, one cell at −0.11). The coordinate-wise chain proposed as a cheaper rescue is refuted separately: a single coordinate's nearest pair matches the pooled nearest pair 2.4 %–6.5 % of the time against a 1.52 % chance rate. clustering::ward_hac is unaffected — it is Euclidean, and its control row is clean at every dimension and size. - Windowed stream queries, and the measurement that says why they are not done by subtraction. CluStream (Aggarwal et al., VLDB 2003) answers "cluster the window [t₀, t₁]" by keeping snapshots and subtracting the one at t₀ from the one at t₁. Cluster-feature additivity makes that exact in real arithmetic, and subtraction is the one CF operation BETULA's stable form does not protect: the entire reason for carrying (n, μ, S) instead of (n, Σx, Σx²) is never to form a difference of two nearly equal large quantities, and an inverse merge forms exactly that, one level up.

window::Moments::checked_subtract implements the inverse merge and conditions it. The two lines

μ_B = (n_AB·μ_AB − n_A·μ_A) / n_B
S_B = S_AB − S_A − (n_A n_B / n_AB)·‖μ_A − μ_B‖²

lose about log₁₀(n_AB/n_B) and log₁₀(S_AB/S_B) decimal digits respectively — and those two ratios are not interchangeable. On a stationary stream S grows with n and they agree. Under drift S_AB picks up the between-window displacement, which has nothing to do with either window's internal spread, and the scatter ratio runs away while the mass ratio does not. Measured on a two-half fixture with a mass ratio of exactly 2.0 — a guard written on point counts sees nothing at all — the recovered scatter is 6.25e-2 against a true 1.0153e-5: a relative error of 6155×, at condition number 4.0e15 and 0.35 surviving digits. That is the finding, and it is why the guard conditions on both ratios a posteriori and returns a SubtractError instead of a number with nothing in it. (A Cholesky downdate restores the definiteness the cancellation destroys — a real and separate problem for the full-covariance feature — but cannot restore digits that were never stored. It treats a symptom of this, not this.)

window::WindowIndex is the constructive answer: micro-clusters stored per frame rather than cumulatively, so a window is a sum of frames and every combination is the stable Chan merge. window_moments and cluster_window answer the CluStream question with no condition number at all — on the same drifting fixture the summed answer is exact to 1e-12 where the subtraction was wrong by 6155×. The trade is explicit and runs the other way from CluStream's: exactness in real arithmetic is given up (a window resolves only to the frame boundary) for an answer that is sound in floating point at every ratio. Capacity is enforced by merging the two oldest adjacent frames, so resolution coarsens with age and never with recency — the pyramidal property, reached by merging instead of by differencing — and compaction preserves total mass exactly, since every step of it is a CF merge.

window::TimeSpan carries timestamp moments under the same (weight, mean, ssd) contract the spatial feature uses, so a windowed summary can say when it is from. Moments is deliberately not a ClusterFeature: the inverse merge is well defined only for the scalar scatter, and letting a diagonal or full feature through it would silently drop the off-diagonal structure a window query never asked to lose.

WindowStream is the estimator on top: partial_fit(X, t) routes points into a CF-tree, seals a frame every frame_width of stream time (snapping forward across a gap so the frame grid stays aligned to the origin rather than to the last arrival), and hands the sealed frame's micro-clusters to the index. close_frame, frame_spans, window_moments and cluster_window are exposed in Python. This is what DenStream structurally cannot do: decay leaves it with only a present, while this retains history and coarsens it with age. - Two benchmark studies the published tables never ran: leaf budget, and insertion order. bench/leaf_budget.py sweeps max_leaves × routing distance × head over three datasets (252 cells, median of seeds 0/1/2) and bench/insertion_order.py measures order sensitivity against a reseeding control (54 cells, 8 runs each). Three findings land in bench/RESULTS.md:

The published digits rows run at max_leaves=4000 against n=1797 — zero compression, confirmed at 1797 leaves, maximum leaf weight 1, mean squared radius 0. They are raw-point clustering behind a betula wrapper and support no claim about summarization. They also understate the library: the curve has an interior optimum, and halving the leaf count raises ward from 0.6428 to 0.6819 and kmeans from 0.4670 to 0.5600. On covtype-20k the ward head holds 0.1412–0.1430 from ×11 to ×202.

The routing distance is a function of compression, mechanically: at ×1.0 the spread across the four distances is exactly 0.0000 on all three datasets, because they build the identical singleton leaf set. Task 27 measured that lever at one budget and found it small — this says the budget it was measured at is the one regime where it provably cannot matter.

And at real compression the input order moves ARI more than the seed does: MNIST at 200 leaves with the k-means head disagrees with itself at pairwise ARI 0.2949 across permutations against 0.7026 across seeds. The ward head's seed arm reads spread 0.0000 / pairwise 1.0000 in all nine of its cells, which is the control that makes the order arm readable. - export_coreset(size=…): the word "coreset" turned into a claim. export_coreset() has always returned the leaf summary and always called it a coreset; nothing checked that it was one. It now carries the two bounds that make the name earned, and they are kept apart on purpose, because the summary and the sample fail in different ways and folding either into the other makes both unfalsifiable.

Summarization, reported in both modes as coreset.offset = Δ = Σᵢ Sᵢ. The summary's cost ĉost(C) = Σᵢ (Sᵢ + nᵢ‖μᵢ − C‖²) is exactly the cost of sending every point of a leaf to the centre nearest that leaf's centroid — the CF identity makes it exact, not an estimate — so it can only over-charge, and boundedly:

0 ≤ ĉost(C) − cost(C) ≤ 4·√(Δ · cost(C)) + 4·Δ         for every C, every k

The proof is one triangle inequality (‖p − c_i‖ ≤ 2‖p − μᵢ‖ + ‖p − c*‖, using only that c_i is nearest to μᵢ) plus Cauchy–Schwarz on the cross term. The relative form is 4√ρ + 4ρ at ρ = Δ/cost(C), and cost(C) ≥ OPT_k bounds it uniformly. Coreset.summary_epsilon(alpha) evaluates it and requires the approximation factor as an argument: reference_cost upper-bounds OPT_k, so Δ/reference_cost under-states ρ and summary_epsilon(1.0) is an optimistic reading, not a certificate. Defaulting that argument would have hidden exactly the assumption the bound stands on.

Sampling, only when size is given. ĉost(C) = Δ + Σᵢ nᵢ‖μᵢ − C‖² and Δ does not depend on C, so the sample only has to be a coreset of the weighted point set {(μᵢ, nᵢ)} — a standard object the published theory covers verbatim — and offset carries the constant rather than losing it. The construction is sensitivity sampling (Feldman & Langberg, STOC 2011) from an α-approximate solution over the leaves, sᵢ = 8·cᵢ/Ctot + 2·(nᵢ/Wⱼ)·(Cⱼ/Ctot) + 4·(nᵢ/Wⱼ), whose total is 10 + 4k — asserted, since a value away from it means the reference solution left a cluster empty. Sensitivity sampling is now known to attain the optimal worst-case coreset size Õ(k·ε⁻²·min(√k, ε⁻²)), matching the STOC 2022 lower bound, and Õ(k/ε²) on Ω(1)-stable instances (arXiv 2405.01339) — which is why it replaced BICO's ring hierarchy in the design.

Coreset.cost(centers) adds the offset for you. That is deliberate: forgetting it understates every cost by the same constant, which is invisible until two coresets are compared.

export_coreset() with no arguments is unchanged and still free — same three arrays, same single O(n_leaves) pass, no k-means. The guarantee numbers it cannot compute are None rather than faked. A size at or above the leaf count returns every leaf exactly instead of a noisy redraw of something already held exactly.

The (k, ε) property is tested the way it is claimed: over candidate solutions the sensitivities never saw, on both sides of the binding. Also asserted — that the summary never under-estimates, that the derived bound holds against the actual points, that the bound is not vacuous on the fixture it is checked on, that the estimator is unbiased over 80 seeds, and that a larger sample is a better one. - The four non-Ward linkages: average, weighted, centroid, median. method="ward" was the only hierarchical head, because the nearest-neighbour chain that drives it is valid only for a reducible linkage — and centroid and median linkage are not reducible, they invert. The four others now run on an Anderberg driver (one global minimum per step over a lazily repaired nearest-neighbour cache; O(m²·d) expected, O(m·d) extra space) and take the same n_clusters, including n_clusters=0 for a Calinski–Harabasz-scored cut. Names follow SciPy's linkage(method=…): average = UPGMA, weighted = WPGMA (McQuitty), centroid = UPGMC, median = WPGMC.

All five linkages are one accumulator. Writing α for the weights a cluster gives its member leaves, each of them is a function of m_A = Σ α_i μ_i and V_A = Σ α_i (S_i/n_i + ‖μ_i − m_A‖²), and both compose exactly under a merge with β = w_a/(w_a + w_b):

m_AB = β m_A + (1−β) m_B
V_AB = β V_A + (1−β) V_B + β(1−β)‖m_A − m_B‖²

The V recurrence is König–Huygens arranged so every term is non-negative — it never forms Σ α‖μ‖² − ‖m‖², which is exactly the cancellation this library exists to avoid, and a test pins it against the CF merge at an offset of 10⁷. The five linkages are then two independent choices: weight the children by mass (UPGMA/UPGMC/Ward) or equally (WPGMA/WPGMC), and measure ‖Δm‖² + V_a + V_b, ‖Δm‖², or Ward's 2·n_a n_b/(n_a+n_b)·‖Δm‖².

On mass weights the accumulator is the merged cluster feature, so three of them are exactly the CF distances the tree already routes by, all on squared distances: UPGMA is D2², UPGMC is D0², Ward is 2·D4². The factor two on Ward is not decoration — it is what puts all five on one scale, where each reduces to the plain squared distance between two points on single-point leaves. Both facts are asserted rather than asserted-about.

weighted and median are the reason the driver is not parameterised over CFDistance: a cluster-feature merge is mass-weighted by construction, so nothing built out of cluster features can represent a cluster whose children were combined equally regardless of their size.

Verification is against implementations that share no algebra with the driver: every linkage reproduces the textbook Lance–Williams recurrence step for step and height for height; centroid and median reproduce SciPy's partition exactly on a one-leaf-per-point tree (average and weighted are excluded there, since SciPy applies them to unsquared distances); and the Ward arm reproduces the existing nearest-neighbour chain at every k. centroid and median are documented as admitting inversions, which is why cuts are taken as a prefix of the agglomeration order rather than by sorting on height. - validity(): three internal indices off the leaf summary, and the caveat each one carries. est.validity() returns calinski_harabasz, davies_bouldin and medoid_silhouette for the fitted partition in O(ℓ·k·d) — no second pass over the data and no O(N²) term, because the sum of squared distances inside a leaf is S_i + n_i‖μ_i − c‖² exactly. On a one-leaf-per-point tree calinski_harabasz reproduces scikit-learn's point-level calinski_harabasz_score to floating-point noise, and a test asserts it.

What is exact and what is a variant is stated rather than glossed. CH is exact. Davies–Bouldin ships the RMS dispersion √(E‖x − c‖²); the classical mean-distance form is not a function of a cluster feature at all, and Jensen only bounds it, so this is a deliberate variant, not an approximation of the original. The medoid silhouette is the index of the summary — a per-leaf ratio weighted by leaf mass, and a ratio of expectations is not the expectation of a ratio. Its medoid is taken in the squared metric, where Σ_{i'} n_{i'}‖μ_i − μ_{i'}‖² = n_j‖μ_i − c_j‖² + const makes the minimiser exactly the leaf nearest the centroid — an O(m) scan where the unsquared medoid needs O(m²).

None of the three can report "there is no structure here." Schubert, Stop using the elbow criterion for k-means (SIGKDD Explorations 25(1), 2023, arXiv 2212.12189), Table 1 has the distance-based indices finding 3–22 clusters in pure noise where BIC correctly finds one. CH is undefined at k = 1, which is that same limitation said out loud. The n_clusters=0 BIC path is unchanged and remains the authority on whether there is one cluster at all. - graph_degree: the density head off its complete graph. method="hdbscan" was quadratic in the leaf count twice over — a full sort per leaf for the core distances, then Prim over the complete mutual-reachability graph — so the one head that most wants a fine summary was the one that could least afford it. graph_degree > 0 (new keyword, default 0 = the exact path, byte-identical to before) replaces both with the two-pass construction of Okkels et al. (Inf. Syst. 142 (2026) 102768, Alg. 4): build a bounded-degree approximate k-NN graph over the leaf means, read the core distances off that graph, take an exact MST of it. Two-pass rather than one-pass because it is the variant whose evaluation reaches cophenetic correlation ≥ 0.9 where one-pass saturates at 0.7, and because its edge count is fixed at O(m · degree) where one-pass lets the expansion queue grow quadratically.

The index is flat — no HNSW layer stack — following Thordsen & Schubert (SISAP 2025), who report that the hierarchy buys little in high dimension, that approximating an RNG/SSG there is futile, and that the capped beam search is the part worth keeping. Three uniformly random out-edges per vertex stand in for the long edges the upper layers would have contributed; they are held outside the degree cap, so a random draw can never evict a true neighbour, and they are what make the diameter logarithmic — without them a capped search on a purely local graph cannot reach a distant region at all.

Median of seeds 0/1/2, one BLAS thread, min_cluster_size = N/100, min_samples = 4N/max_leaves, timing the head only (the tree build is identical and subtracted): on MNIST (70 k × 784) at max_leaves=8000 the head goes 52.0 s → 0.45 s at degree 32, a 116× cut with the ARI unchanged at 0.0523 and the same ten clusters — the approximation is free there. On 100 k 2-D blobs at max_leaves=32000, 36.6 s → 0.98 s (37×) for ARI 0.5674 → 0.5608; on covtype (581 k, 54-D) at max_leaves=8000, 4.61 s → 0.45 s for ARI 0.0490 → 0.0457. The trade is monotone in the degree and it is the degree, not the graph, that costs the quality — degree 8 is not enough at d = 784 (ARI 0.0190 with 4 clusters, swinging across [0.0068, 0.0600] between seeds), 16 is the smallest measured to be lossless there.

graph_degree is a floor, not a ceiling: core distances read off a graph saturate at the farthest neighbour, so a degree below what min_samples needs underestimates every core distance with no bound on the error. The head raises the request to min_samples / mean leaf mass whenever that is larger — the paper's Ω(minPts) requirement translated into the currency a weighted summary actually counts in. - method="mppca": a mixture of probabilistic PCA, and the compression trade it exposes. Each component covariance is Σ_c = W_c W_cᵀ + σ_c² I with W_c of rank rank (new keyword, default 2, clamped to dim - 1; 0 gives a spherical mixture). Both identities it rests on — the Woodbury inverse (1/σ²)[I − W M⁻¹ Wᵀ] and |Σ| = σ^(2(d−q))|M|, M = σ²I_q + WᵀW — were verified exactly in Maxima at three shapes before any of it was written, and nothing in the head or in point scoring ever forms a d×d matrix. That is what makes it run at d = 784 and max_leaves=2000, where gmm-full's per-leaf dense scatters need ~38 GB across its parallel restarts and the benchmark table has carried a — for it since 0.1.0.

At full resolution it beats the diagonal and the full head. On digits (1797×64, feature="fd", max_leaves=2000 ⇒ one leaf per point, median of seeds 0/1/2) rank=5 scores ARI 0.600 against gmm-full's 0.575 and the diagonal head's 0.461 — a rank-5 subspace where full covariance spends 2016 parameters. On six 5-dimensional subspaces sharing one centre in 100-D, where centroid distance carries no information at all, ARI peaks exactly at the true rank: 0.385 / 0.654 / 0.998 / 0.823 / 0.727 for rank 2 / 3 / 5 / 10 / 20, against the diagonal head's 0.166 (median of seeds 0/1/2; the band at the peak is [0.9976, 0.9984]). Pull the same six subspaces apart so the centroids separate them and every head scores 1.0000 at every one of those ranks — the extra parameters cost nothing where they buy nothing.

And it loses on MNIST, for a reason that is about the summary rather than the head. The expected-log E-step folds each leaf's own scatter into the component covariance as −½ tr(Σ_c⁻¹ Σ_i), and that within-leaf term carries almost none of the between-cluster orientation — so compression costs a head in proportion to how much orientation it models. Sweeping the leaf budget on digits shows the ordering invert cleanly: at one leaf per point mppca(5) 0.600 > gmm-full 0.575 > gmm 0.461; at 6:1 compression gmm 0.493 > mppca(5) 0.406 > gmm-full 0.273; at 16:1 0.235 > 0.168 > 0.099. MNIST-20k at max_leaves=2000 is 10.6:1, and there mppca scores 0.159 / 0.069 / 0.024 at rank 2 / 5 / 10 against gmm's 0.274. The head is opt-in, the default is untouched, and docs/USAGE.md states the loss and the mechanism.

predict / predict_proba work as for any mixture head, and n_clusters=0 selects the component count by BIC (the Stiefel rotation of W_c is unidentifiable, so q(q−1)/2 of the loadings are not counted as free parameters). - absorb="subspace": the χ² gate read on the leaf's own low-rank basis. Only feature="fd" carries a basis, so every other feature model falls back to chi2 and the option changes nothing unless you asked for the Frequent-Directions sketch. Same chi2_p and chi2_scale, same units. Under the same Normal-Inverse-Gamma prior the effective covariance is Σ_eff = a·FᵀF + b·I, and Woodbury reduces its inverse to one ℓ×ℓ solve — no d×d matrix is ever formed. Gram-Schmidt on the sketch rows would be cheaper and would be wrong: FD rows are orthogonal only immediately after a reduce(), and leaves are queried in between.

It works where orientation is the only signal, and loses on MNIST. On six rank-5 subspaces sharing one centre — centroid distance carrying no information whatsoever — leaf purity goes 0.820 → 0.938 (median of seeds 0/1/2, max_leaves=2000, chi2_scale=0.01, seed ranges disjoint), and on separated blobs it reaches the same ARI 1.0000 with 6 leaves against 99, a 16× smaller summary at identical quality. On MNIST-20k (784-D, StandardScaler, k=10) it scores ARI 0.250–0.260 against chi2's 0.274–0.291 at every scale tried, with more leaves and ~20 % more time. Opt-in, default untouched, and docs/USAGE.md states the loss.

The bound worth knowing is not the gate's. On that concentric fixture both gates reach purity 0.82–0.96 and ARI ≈ 0.05: every head here assigns by centroid, and there every centroid coincides. A better-oriented summary currently has nothing that can read it, which is why this is not a candidate default. - projection="svd": a CF-weighted PCA of the leaf summary, and the one-call text pipeline. Clustering TF-IDF in its own geometry does not work — the sparse path scores ARI 0.003 on 20-newsgroups — and the standard fix is to reduce first. This does the reduction inside the same call, factorizing the M ≈ 10³ leaf centroids rather than the N documents: on 20-newsgroups (18 846 × 2 000, k=20, rank 50, median of seeds 0/1/2) fit_predict_sparse(..., projection="svd", max_leaves=256) scores 0.130 in 0.30 s, against TruncatedSVD(50) + KMeans at 0.143 in 0.54 s; at max_leaves=2048 it reaches 0.152 in 5.4 s. Available on fit_predict, the Betula estimator and fit_predict_sparse, and unlike the NMF projection it accepts signed data.

The projection is only half of it — the other half is that a PCA is a linear map. A raw row encodes as (x − x̄)Vᵀ, so a projected fit keeps the head's own point rule instead of falling back to "answer with the row's leaf's label". Measured on the same basis and the same rows, that distinction alone is worth 0.062 ARI (0.097 → 0.159). The NMF projection cannot be given the same treatment: its code is the solution of a per-row nonnegative least squares, not a matrix product, so it keeps the microcluster route and the docs now say why.

Two measured facts decide whether it works for a given corpus, both in docs/USAGE.md. Use a cosine head: method="kmeans" on the same codes scores 0.014 against spherical-kmeans's 0.152, because the leading principal direction of a TF-IDF corpus is document length. And the leaf budget is the cost, not the projection: sweeping rank 1 → 100 moves the total by 1.2 s, while max_leaves 256 → 2048 moves it from 0.30 s to 5.4 s, the sparse summarizer being quadratic in the micro-cluster count.

The basis is not a compromise for being built from a summary. Labelling raw rows in it scores 0.159 against 0.143 for TruncatedSVD's own basis on the same rows — under the spherical cluster feature the discarded within-leaf scatter is isotropic, so it shifts eigenvalues and leaves the directions alone. bench/_worker.py's betula-svd row now runs betula's own reduction; it previously borrowed scikit-learn's TruncatedSVD and so measured scikit-learn's reducer. - refine=n: BIRCH's Phase 4, which this crate never had. n Lloyd sweeps over the raw rows, warm-started from the Phase-3 centres. Off by default, centroid heads only (kmeans, spherical-kmeans — a mixture assigns by maximum posterior and Ward/Spectral/Leiden by microcluster, so sweeping centres would substitute a partition they never fit), and in-memory fit / fit_predict only: partial_fit keeps a tree rather than the data, and the CSR path would have to densify the matrix it exists to avoid. Label-changing where enabled — which is why it is opt-in rather than a new default.

It moves the objective exactly where the theory says — a summary coarse relative to the data — and on this benchmark that is not the same as moving the answer. On MNIST (first 20 000 rows, 784-D, StandardScaler, k=10, spherical CF, max_leaves=4000, median of seeds 0/1/2) twenty sweeps take the k-means objective 11 750 563 → 11 710 630 for 4.2 s → 6.4 s, and take ARI 0.315 → 0.309 with it. Five sweeps: 11 720 402 and 0.311. The descent is monotone in the objective and monotone downwards in ARI.

The acceptance target was not met, and the misses are the interesting part. Against sklearn.cluster.KMeans(n_init=10) (ARI 0.324, objective 11 671 351, 19.3 s) refinement lands at 3.0× the speed but 0.34 % above the objective: one warm start from a lossy summary settles in a worse basin than the best of ten k-means++ restarts, and more sweeps cannot leave a basin. digits and covtype gain nothing at all, structurally — digits at max_leaves=4000 realises 1 797 leaves for 1 797 rows, so Phase 3 already is exact k-means on the raw points and Phase 4 begins at its fixed point, while covtype moves 0.1993 → 0.1998. MNIST at max_leaves=16000 behaves like digits: 0.3237 at objective 11 671 813 — scikit-learn's own answer in 15.8 s against its 19.3 s — and twenty sweeps then move the ARI by 0.0001.

The measurement also sharpens the caveat the feature ships with: on this benchmark a lower k-means objective goes with a worse partition, twice. covtype — n_init=10 objective 827 314 / ARI 0.174, n_init=1 objective 832 081 / ARI 0.277. digits — n_init=10 69 405 / 0.468, n_init=1 69 749 / 0.559. refine optimizes the objective faithfully; docs/USAGE.md says plainly that this is not the same as optimizing the answer. - The full BIRCH absorption grid, and D3 implemented for the first time. absorb accepted only euclidean and chi2; it now takes manhattan (D1), average (D2), diameter (D3), ward (D4) and radius (R) as well. D3 — the mean squared distance within the cell that results from the merge — did not exist in src/distance.rs at all, so a third of the criteria a BIRCH paper refers to were unreachable.

D3 closes over (n, μ, S) because the double sum telescopes: Σᵢⱼ‖xᵢ−xⱼ‖² = 2·n·S, leaving D3²(A,B) = 2·(S_A + S_B + n_A·n_B/(n_A+n_B)·‖Δμ‖²)/(n_A+n_B−1) — D4's merge term and the two scatters over one fewer than the merged mass. Checked against brute-force enumeration of the underlying points, both in a Rust fixture and over random point sets to 5.7e-14.

The default does not move, and docs/USAGE.md now says why rather than leaving it implicit: the thesis tunes absorption for minimum variance and prefers D4/D2, while this crate chose mass-invariance, and the two objectives genuinely conflict — the variance-minimising criteria inherit the BIRCH size-imbalance bug that chi2 exists to fix. threshold is read in each criterion's own units (L1 for manhattan, squared for the rest, a χ² quantile for chi2), so it does not transfer between them; the table in docs/USAGE.md states the units per criterion.

The one-shot path had carried a verbatim copy of the estimator's gate resolution, which would have meant adding every criterion twice; it now calls the same resolve_gate, and the accepted-value list exists once so the parser and the error messages cannot drift apart. - A UserWarning when the CF-tree compressed nothing: one leaf per point. With threshold = 0 a point is absorbed only by an entry it equals exactly, so a max_leaves that never binds leaves the summary at n micro-clusters. The answer is not wrong — but the tree is descended, split and rebuilt for every point and Phase 3 then runs on the raw rows anyway. Measured single-threaded on 64-dimensional blobs, that costs 3.8× the fit time at n = 8 000 and 14× at n = 40 000 against the same call with a binding max_leaves = 2000, and 36–46× more than clustering the raw rows directly. The default max_leaves = 2000 means only a caller who raised it explicitly can land here, which is exactly why it was silent.

The check reads the realised leaf count rather than comparing n against max_leaves, matching the existing budget warning's own rule: leaves == n is the thing that is actually true, and only it says no two points ever shared a micro-cluster.

This closes task #40's third item by refuting it. The proposed fix was a leaf-per-point bypass that skips the routing descent in this regime; it cannot be an exact equivalence. At threshold = 0 an exact duplicate is absorbed — but only when it descends to the leaf holding its twin, and which leaf a point reaches depends on the tree's state at the moment it was inserted. A hash-based bypass would merge duplicates the real tree keeps apart, so it would be a relabelling sold as an optimisation. The measurable cost is published and made visible instead.

• A UserWarning when the CF summary is too coarse to carry n_clusters. Asking for k clusters from fewer than 2k leaves silently produced a worse partition with no signal; the warning names the realised leaf count, k, their ratio and the current max_leaves, and points at the three parameters that fix it. Fires from every entry point that partitions a summary into a caller-supplied k — the fit_predict function, the Betula estimator (dense and CSR), fit_predict_sparse and KPrototypes — and stays silent for auto-k and for the heads that discover their own count (hdbscan, scale-space, leiden).

The threshold is measured, not assumed: over three seeds on the ward head, well-separated data is already at its plateau ARI at ≈2 leaves per cluster and loses 29 % (k=50) / 55 % (k=200) of it at ≈1, while digits and covtype score 0.000 and 0.003 at ≈1. Lang's thesis reports the same floor from the other end (Sec. 5.5.4). It is a floor, not a recommendation — covtype peaks at ≈8 leaves per cluster and declines after — and the check reads the realised leaf count, since the tree routinely settles below its cap and at n < max_leaves the cap never binds at all.

• A warning when a Gaussian head is fed a feature that cannot carry its covariance. method="gmm" and method="gmm-full" fit a per-dimension covariance, but Spherical::variance(_d) ignores its argument — it returns ssd / (w · dim), one isotropic number for every dimension, because a spherical cluster feature holds a scalar scatter and cannot hold more. The M-step adds that number to all dim component variances (and gmm-full inherits it through cov_dense's diagonal default). With one leaf per point ssd = 0 and nothing is added; under compression each component is inflated equally in every dimension by however much leaf scatter it happens to cover, so a dimension with genuinely near-zero variance is lifted to the isotropic average.

The fit survives that. The labelling does not, because the maximum-posterior argmax is dominated by ln|Σ_c| = Σ_d ln σ²_cd while a nearest-centroid rule ignores Σ. Measured on digits, medians of seeds 0/1/2: at 898 leaves (×2.0) the head returns ARI 0.0088, while labelling the same points off the same fitted model by nearest centre returns 0.5288 and by their leaf 0.5296. Ten non-empty components, largest holding 30% — there is no degenerate component and no merged blob.

The falsifiable prediction — features carrying per-dimension scatter must not show it — holds: fd and full read 0.3840 and 0.4403 in that same cell, and all three agree to the digit at ×1.0, where there is no scatter to add. gmm-full on the spherical feature collapses identically (0.0096 at 1200 leaves) and never does on feature="full".

So the pairing is now a warning naming the measured cost and the fix (feature="full", or "fd" in high dimension). This closes task #89 with a mechanism rather than a workaround; heads other than these two read the same isotropic variance(d) but were not measured, so the warning does not claim them. - scripts/check_mutants_baseline.py — the staleness half of the mutation ratchet, in a second rather than three hours. mutants-baseline.txt records accepted surviving mutants by their exact path:line:col: description, so any edit that shifts a line silently invalidates every entry below it in that file. Nothing checked that: the weekly mutants.yml run notices only as a "trim them" note inside a step summary, and it had been red since 2026-08-19 for an unrelated reason. Checked by hand on 2026-09-07, 147 of the 315 entries no longer matched anything — 121 the same mutation at a new line, 26 gone entirely, worst in stream.rs (33), scalespace.rs (28) and hdbscan.rs (21, rewritten for the parallel build).

cargo mutants --list builds and runs nothing — about a second on this crate — so this is answerable on every push. The script separates the two cases that need different work: a moved entry keeps its recorded argument and needs re-anchoring, and 42 of the 121 are unambiguous enough for the script to propose the new line outright; a vanished entry describes code that is gone and takes its justification with it. It is not yet a CI job.

It also surfaced a second problem it was not written for: the crate now lists 8 927 mutants against the 4 312 that .cargo/mutants.toml sizes the 96-shard matrix on, so the weekly run cannot finish as configured even once it is green again. (The tagged v0.8.0 text reached that conclusion through wrong arithmetic — it multiplied a stale 101 worker-seconds per mutant and forgot -j4. The measured basis is in [Unreleased]; the conclusion survives it.)

Changed

• The sparse path stores micro-cluster coordinates feature-major, and the 20-newsgroups fit is 3.4× faster with every label bit-identical. Both sparse passes — the leader summarisation and the row labelling — held each micro-cluster's coordinates as its own dense Vec<f64> of length n_features. Scoring one row against L of them therefore pulled L · nnz scattered cache lines out of L · n_features · 8 bytes of state, none of it reused before eviction: 32 MB at L = 2048 on the 20news fixture. The symptom was an elbow in the leaf-budget curve exactly where that working set leaves L3 — doubling max_leaves from 512 to 1024 cost 6.3×, not 2×.

Transposed to [c · L + i], a row instead walks nnz contiguous runs of L doubles, which is both cache-resident and an elementwise axpy the compiler may vectorize — unlike the per-cluster dot product, which is a float reduction and so may not be reassociated. The arithmetic is deliberately untouched: for a fixed cluster the row's terms are still accumulated in row order, so every distance, absorption decision and scatter update is bit-for-bit what the old layout produced. That is measured, not asserted — all twenty label digests across five budgets, three seeds and both pipelines are identical between the two builds.

Measured single-threaded on 18 846 × 2 000 TF-IDF: betula-sparse 6.57 s → 1.94 s, betula-svd 5.12 s → 3.31 s, and at max_leaves = 2048 on the raw kmeans path 12.76 s → 2.19 s (5.8×). The super-linear elbow is gone. Peak RSS is unchanged (332 MB): the transposed buffer is bounded by min(max_leaves, n_rows) · n_features, the same worst case as one vector per cluster, and it widens in place rather than through a second live copy.

A perf profile is what made this complete rather than half-done. After fixing the summarisation pass named by the task the fit was only 1.5× faster, and the profile put 65% of the remaining samples in the row-labelling pass against 6% in summarisation. The task named one half; the profile named the other.

Breaking (Rust API): sparse::nearest_sparse(&means, &musq, …) is replaced by sparse::SparseCentroids::from_features(&micros).nearest(…). The free function took the centroids and their cached ‖μ‖² as two independent slices that had to agree, and nothing could enforce it; the type derives one from the other. The Python API is unchanged. - The FAISS baseline is now measured in two configurations, and the "quality-per-second" hedge is withdrawn. bench/external_baselines.py timed faiss.Kmeans at its defaults, which is not a like-for-like fit: max_points_per_centroid defaults to 256, so at k = 8 FAISS trains on 2 048 of the 200 000 rows and assigns the rest, and it seeds from a random subset rather than k-means++, at nredo = 1. Two of the three things that made the row fast are quality-costing shortcuts.

A new faiss-kmeans-matched row hands it every row and ten restarts — the cheapest setting whose median ARI reaches betula's. It costs 1.19 s at 200 000 (6.3× betula) and 5.38 s at 1 000 000 (7.7× betula), and at 1 M still returns 0.835. The ceiling is initialisation, and scikit-learn is the control that proves it: on the same fixture with every row used, init="random" medians 0.692 at n_init=1 and 0.835 at n_init=10, while init="k-means++" medians 1.000 at n_init=1 on every seed. ARI 0.835 is a persistent local optimum — two blobs merged, one split — that neither more data nor more iterations escapes: ARI is flat from niter = 25 to niter = 300.

So task #73's premise was wrong. There was no throughput gap to close; betula's edge is k-means++ on the exact CF potential over a summary small enough that restarts are free. The board totals are unchanged (8/54/6 quality, 32/1/5 speed, 30/4/0 memory) and the two default-configuration FAISS speed losses stand as measured. - method="ward" with n_clusters=0 selects k by the variance ratio, not by the elbow. The auto-k path cut the dendrogram at the largest relative jump in merge height, which is the elbow criterion wearing a dendrogram, and Schubert (SIGKDD Explorations 25(1), 2023) is a direct argument against it. It now scores every horizontal cut in the sweep by Calinski–Harabasz — the first alternative that paper names, exact on cluster features — and keeps the best. On two far groups of two nearby subclusters each, the tallest relative jump is the one that joins the far groups, so the old rule reported k = 2 on every seed where the new one reports the true 4; that fixture is now a test, and it fails if the elbow is put back. Scoring 19 cuts is O(k_max·m·d), nothing beside the O(m²) dendrogram that had to be built anyway. k = 1 is still outside what CH can express — n_clusters=0 on a mixture head remains the way to ask whether there is structure at all. - spherical-kmeans skips the dot products it can prove will not change a label. Hamerly's two bounds, restated on similarities rather than distances (Schubert, Lang & Feher, arXiv 2107.04074): keep a lower bound low[i] on the cosine to the assigned centre and an upper bound high[i] on the best of the others, loosen both by how far the centres turned, and skip the leaf entirely while low[i] ≥ high[i]. Labels, centres and iteration count are unchanged — byte-identical SHA-256 label digests in all nine measured cells, and the Rust suite checks the bounded loop against an independent re-derivation on both fixtures.

Restating the bound on the sphere is not a term-by-term transcription, and the place it breaks is worth naming: cos(θ − φ) is an upper bound on the reachable similarity only while θ ≥ φ. A centre that turns by more than the angle it started at sweeps through the point and reaches similarity 1, so the aggregation Hamerly gets for free — nearest of the other centres paired with the largest movement among them — is unsound here, because d − δ is monotone in d and cos(θ − φ) is not monotone in θ. The symmetric failure sits on the other side: past θ + φ = π the cosine turns back upward and cos(θ + φ) stops being a lower bound. Both are guarded, both guards have a test that fails when the guard is removed.

Median of seeds 0/1/2, one BLAS thread, timing the head alone (max_iter=100 minus max_iter=1, which never reaches the skip test, so the tree build and the first assignment pass cancel):

data	leaves	k	head, no skip	head, bounded	speed-up
100 k × 20-D blobs	32000	10	0.28 s	0.09 s	3.04×
100 k × 20-D blobs	32000	50	3.55 s	1.55 s	2.29×
100 k × 20-D blobs	32000	200	12.01 s	7.34 s	1.64×
MNIST 70 k × 784-D	8000	10	3.41 s	2.61 s	1.31×
MNIST 70 k × 784-D	8000	50	13.82 s	11.27 s	1.23×
MNIST 70 k × 784-D	8000	200	82.76 s	68.53 s	1.21×

The bound is a dimension effect, not a k effect, and the table says so twice over: it is worth 3× at d = 20 and 1.2–1.3× at d = 784, and within each dataset the gain falls as k grows. In high dimension the cosines from one leaf to every centre concentrate, low ≥ high almost never holds, and what is left is the bookkeeping. This is the same concentration of measure that makes the Euclidean path collapse on MNIST — it is not a defect of this implementation, and no amount of tightening will make Hamerly pay there. method="vmf" shows no gain outside measurement noise (1.05× against a base-measurement spread of the same size), which is the expected answer: its E-step is soft, every leaf needs all k responsibilities, and only the one-restart spherical k-means warm start is accelerated at all. - The CF k-means++ sampling weight now carries the leaf's own scatter. Seeding drew each candidate centre proportional to n_i·D²_i — the point-level k-means++ weight, with the leaf treated as a point at its mean. The exact CF-adapted potential is Lang's Eq. 5.4, S_i + n_i·D²_i, and the missing S_i is not a rounding detail: a wide leaf sitting on an already-chosen centre has D²_i = 0 and so was unsamplable however much scatter it held. The greedy candidate score inside the trial loop is unchanged and correct as it was — it differs from the exact potential by Σ_i S_i, which is constant in the candidate. cop_kmeans gets the same weight, over a chunklet scatter assembled by König-Huygens (Σ S_i + Σ n_i‖μ_i − c‖²), so with no must-links it seeds identically to kmeans.

This is on correctness grounds and claims no quality gain — it relabels output and the benchmark moved in both directions. It reaches kmeans, gmm, gmm-full, spectral, xmeans and cop_kmeans; ward, leiden, hdbscan, vmf and spherical-kmeans (which has its own angular seeder) are untouched, and every scikit-learn row reproduced to the last digit. It is a no-op wherever the summary is lossless: at max_leaves ≥ N every leaf holds one point, S_i is 0, and the weight is bit-identical to the old one — which is why the whole synthetic quality table and every digits row are unchanged, and it is a useful self-check that a benchmark budget actually compresses.

Where it does bite, over three seeds (bench/RESULTS.md re-measured whole): covtype spectral 0.037 → 0.100 and MNIST spectral 0.155 → 0.203; MNIST kmeans 0.302 → 0.307; covtype kmeans 0.088 → 0.074, gmm 0.088 → 0.076. At 16 000 leaves covtype gmm goes 0.094 → 0.104 and kmeans 0.078 → 0.067. The covtype GMM result the README quoted at the default budget was a win by 0.008 and is now a loss by 0.004 — both sit inside two nearly coincident seed ranges (0.055–0.096 against 0.055–0.102), so the honest reading is that it was a tie before and is a tie now, and the page says so. The claim that survives is the one at adequate leaf resolution, which got stronger.

It also inverted the refine= MNIST result above: with the new seeding Phase 3 already lands at ARI 0.315 rather than 0.275, and twenty Lloyd sweeps now lower the objective while lowering the ARI. That entry was re-measured rather than left standing.

A prior sweep of ELKI's three CFInitWeight options at 18×/6×/2× compression over five seeds found no consistent winner and every difference inside the seed spread (local/scratch/elki/init_weight_ab.py); Lang's thesis separates them only at RMSD ≈ 0.005 with ±0.004 confidence intervals. This change is made because the formula is the right one, not because it scores better. - min_samples now counts the microcluster itself in method="hdbscan", so min_samples=1 leaves every core distance at 0 and HDBSCAN* degenerates to single linkage. The convention is a genuine split in the field — Campello, Moulavi & Sander's Def. 3.1, sklearn.cluster.HDBSCAN and ELKI (whose parameter reads "including this point") all include the object, while scikit-learn-contrib/hdbscan excludes it — and it was previously stated nowhere: not in the Rust docs, the .pyi or docs/USAGE.md. The shipped behaviour was the exclusive one, so min_samples=5 acted like the standard min_samples=6, and the published benchmark compared our effective 11 against sklearn-hdbscan's 10.

Aligning with the majority — and with the library this project mirrors and benchmarks against — makes the comparison like-for-like. This relabels method="hdbscan" output: to keep the old behaviour, add one to min_samples. The convention is now pinned by a test on the core distances themselves rather than inferred from a downstream partition.

The four quality tables were re-measured over three seeds after the change. Cell by cell against the previous run, only betula-hdbscan rows moved — every other method reproduced to the last digit on every metric — so the deltas are attributable to this and nothing else. digits goes from a loss to a win (ARI 0.146 → 0.164 against sklearn-hdbscan's 0.149, both zero-spread, with the noise fraction down from 0.620 to 0.580); blobs 0.142 → 0.154 and varied 0.519 → 0.536; aniso 0.569 → 0.568 and covtype +0.0001 are ties; MNIST stays all-noise.

• The three distance kernels get a hand-written AVX2 + FMA path: 1.38x on a 20-dimensional kmeans fit, 1.59x on a 784-dimensional gmm fit, byte-identical labels. src/kernels.rs said its reductions were "plain inlinable loops the compiler vectorizes at each call site". perf annotate on the symbol holding 31.5% of a 1 000 000 x 20 kmeans profile returns only subsd, mulsd and addsd -- not one packed instruction -- so the claim was wrong, and the doc is corrected along with the code.

It is wrong for a reason no build flag can fix. a.iter().zip(b).map(...).sum() is Iterator::sum, a strictly-ordered left fold; IEEE addition is not associative; LLVM may therefore not reassociate the reduction, and without reassociation there is nothing legal to pack.

Measured on this machine against the shipped scalar fold: four independent accumulators at the crate's baseline ISA give 0.93x at d = 20 -- a loss -- 1.13x at 64 and 1.46x at 784, while explicit AVX2 + FMA gives 2.35x / 3.89x / 2.43x. Breaking the dependency chain is not the win; the register width is, and at a baseline x86-64 target the only way to reach it is to write the intrinsics and detect the feature at run time.

So sq_euclidean, dot and manhattan now dispatch to hand-written AVX2 kernels for f64 and f32 when the CPU reports avx2 and fma, the vector is at least 16 wide, and fall back to the identical scalar fold otherwise. Non-x86_64 targets and pre-AVX2 x86_64 are unaffected -- the fallback is the old code, and the choice is made at run time, so one wheel still serves every machine.

The width gate is the part worth reading, because without it this was a regression. A #[target_feature] function cannot be inlined into a caller without the same features, so taking the packed path swaps an inlined loop for a real call plus a horizontal sum. At d = 20 that is noise; at d = 2 the scalar body is two multiply-adds and the fixed cost is the whole function. Ungated, kmeans on 500 000 x 2 measured 0.65x -- 1.55 times slower -- and 0.89x at d = 8. Two-dimensional input is not hypothetical: the crate's own published scaling table is measured on d = 2 blobs. The gate plus #[inline(never)] on the dispatch helper (both were needed; with the helper inlined, d = 8 came out 1.26x slower than the ungated version it was meant to fix) leaves: 0.93x at d = 2, 0.92x at 8, 1.08x at 16, 1.24x at 20, 1.57x on 20 000 x 784 gmm, and 1.18x / 1.58x on float32 at d = 16 / 24. An unexplained ~8% residual below d = 16 is recorded in ADR 003 rather than argued away; a three-build control puts build-to-build noise at +/-1%.

On float32 the label digest does change at d >= 24 -- but the ARI between the two label vectors is exactly 1.000000. The partition is identical; only the cluster numbering permutes, because k-means++ draws its seeds in a different order when a sampling weight moves by an f32 ulp. f64 digests are identical on every shape measured. A digest comparison alone would have called this a behaviour change and an ARI alone would have hidden that the ids move; both numbers are reported because they mean different things.

This is the crate's first unsafe, about 130 lines in one leaf module, and it is a decision rather than a detail: ADR 003 records the alternatives, including the safe four-accumulator version that the numbers rejected. Two preconditions are upheld one call frame above and both are tested -- the feature detection, and n = a.len().min(b.len()), which matters because CFTree::insert accepts a point longer than the tree's dimension and an intrinsic loop bounded by a.len() would read out of bounds rather than merely return a wrong number.

The acceptance gate was that nothing observable changes. Eight partial sums and an FMA round differently from one serial chain, so an argmin decided by less than an ulp could in principle flip, and identical labels are therefore a measurement and not a theorem. What was measured: SHA-256 label digests identical on both profiled shapes across an A-B-A-B run; the whole Rust suite, SciPy and ELKI cross-checks included; and a re-run of the entire quality benchmark, 75 of whose 78 rows reproduce the committed per-seed table cell for cell. The three that moved are all betula-hdbscan, and a direct A/B on exactly those cells returns identical label digests from both builds -- they predate this change and belong to 835d05f, which made min_samples and min_cluster_size count points rather than leaves and relabels the summary route by design.

• method="gmm" hoists the log-variance term out of the leaf loop: 1.35x on the profiled shape, byte-identical labels. A profile came first, as the task required. On 20 000 x 784 with method="gmm" and max_leaves=2000 (perf record -F 997 --call-graph fp, 9 014 samples, RAYON_NUM_THREADS=1), gmm_diagonal_once was 34.9% of the run and __ieee754_log_fma alone was 14.6% -- one library call, a seventh of the whole fit.

The cause is a loop-invariant in the innermost E-step statement: 1/2 log(2 pi sigma^2_cd) and log w_c depend on the component c and the coordinate d, never on the leaf i, but they sat inside the for i in 0..m loop that runs around them. At m = 1833, k = 10, d = 784 that is 14.4 million ln calls per iteration instead of 7 850. Hoisted to the top of the iteration the cost becomes k*(d + 1).

The accumulation order is untouched -- the same addends are summed in the same sequence, only the transcendentals are precomputed -- so this is not an approximation and the acceptance gate was that it changes nothing. Measured A-B-A-B on the same build with only the hoist differing: 2.847s / 2.098s / 2.781s / 2.072s, 1.35x, with the SHA-256 label digest e887117b6aaa25c2 identical in all four runs.

The same sweep over gmm-full, gmm-toeplitz and mppca found the identical pattern -- but there only log w_c is left inside the leaf loop, k calls per leaf against an O(d^2) Mahalanobis term, so hoisting it would buy nothing measurable and those loops are unchanged. Recorded here so the next reader does not re-derive it.

• Every published table re-measured on the current build; the hdbscan rows moved both ways. The committed CSVs predated 835d05f, the commit that made min_samples and min_cluster_size count points rather than leaves and whose own message says the summary route relabels. Nothing else in the record had caught up. Re-run over seeds 0/1/2: only betula-hdbscan quality cells moved, every other method reproducing cell for cell, which attributes the deltas to that commit and nothing else.

It is a trade, not a win, and both halves are published. mnist at n = 20 000 goes from all-noise (ARI 0.000, zero clusters) to 0.117 with 49 clusters, and covtype from 21 clusters to 43 at a flat ARI -- the catastrophic case is fixed. But varied drops 0.536 → 0.479 with the cluster count going 2 → 73, and blobs 0.154 → 0.142: counting in points makes the size gate roughly N/max_leaves times weaker at the benchmark's own settings, so the head now over-splits where it used to under-split. The three-seed range on aniso is 0.016–0.993, which is the honest summary of that row -- a regime without a stable answer rather than six numbers.

Speed moved too, and not uniformly. Streaming 20-D is 1.29× faster and the sparse pipeline 1.10–1.42×; the d = 2 scaling fixture is 0.92×, the AVX2 gate's residual cost on the one shape in the suite where packed kernels cannot help; betula-hdbscan is 0.79–0.84×, the price of the same units fix doing more work per core distance. Full covtype now clusters in 1.06 s against scikit-learn's 4.92 s.

The scoreboard ratchet fired and was accepted deliberately rather than overwritten: one genuine demotion (betula-nmf speed, win → tie, because scikit-learn's NMF got faster by more than ours did) and six vs-best pairings that vanished rather than lost, because the non-betula champion changed identity at 1 M. bench/RESULTS.md names each one.

Fixed

• fit_predict_sparse(method="vmf") returned a single cluster. Labels change on every sparse fit with a directional head. A von Mises-Fisher density lives on the unit sphere, and the estimator knows it: Betula._build forces normalize=True for vmf / spherical-kmeans. But fit_predict_sparse is a module function that never passes through _build, so it handed the head raw TF-IDF rows. leaf_means states the precondition it needs — on L2-normalized data a leaf mean is a resultant of length R̄_i ≤ 1 — and violating it makes ‖μ_i‖ a per-leaf multiplier on κ_c. That is a positive constant only while κ is shared, which is exactly why the same omission in spherical-kmeans is invisible (its argmax is a single similarity, and the norm cancels); with a fitted κ per component it reweights the components against each other, and the Mixture the head hands back — which normalizes, as vMF requires — ends up describing a different model from the labels the head returned. Measured on a 6 000 × 4 000 block-topic corpus at max_leaves=2048, converged after two EM iterations: κ = [3558, 3514, 3552, 1402], the head's own labels [502, 499, 503, 544], and the mixture reassigning all 2048 leaders to the one broad component — 1504 of 2048 leaves disagreeing with the fit they came from. Since 0.6.0 labels every mixture head by maximum posterior, that surfaced as one cluster for the whole corpus: ARI 0.000 against the dense tree's 1.000. The sparse entry point now L2-normalizes its rows for the directional heads, in O(nnz) — scaling the stored values is the whole operation, an implicit zero stays zero. The same fixture reads 0.995; spherical-kmeans is unchanged at 1.000 on both sides, as the cancellation argument predicts. Two regression tests, each failing on revert.

• The sparse path labelled the posterior and agglomerative heads by nearest micro-cluster, which answers a different question. Labels change for gmm, vmf, ward and the other agglomerative heads on sparse input. Two rules replace it. gmm and vmf now label by maximum posterior, the rule the dense path has used since 0.6.0: a diagonal Gaussian's quadratic form splits as Σ_j (x_j − μ_cj)²/σ²_cj = Σ_{j: x_j ≠ 0}(x_j² − 2x_jμ_cj)/σ²_cj + Σ_j μ²_cj/σ²_cj, and the second term is one number per component built once — so the density costs O(nnz·k), not the O(d) this path exists to avoid. A vMF density already reads only the non-zeros. The full-covariance, Toeplitz and probabilistic-PCA kernels have no such split (a Cholesky solve and a Woodbury projection both read every coordinate) and are refused rather than approximated. Everything else — the agglomerative heads, which have no point model at all — now keeps the label of the micro-cluster the summarisation put the row in, which summarize_sparse already computed and threw away: it agrees with the summary by construction, it is the sparse counterpart of the dense path's point-to-leaf route, and it is free where re-deriving the nearest micro-cluster cost O(nnz·L).

On the 6 000 × 4 000 block-topic corpus at max_leaves=2048, median of seeds 0/1/2: gmm 0.026 → 0.110, ward 0.022 → 0.068. Both are still bad, and the reason is the summary, not the rule. The leader pass spends its budget at row 2048 and force-assigns the remaining 3 952 into 544 leaders, so three quarters of the mass sits in mixed micro-clusters; a head that models each one separately inherits that, while kmeans (0.976) and spherical-kmeans (1.000) survive because pooling thousands of leaders into k centroids averages it away. The dense tree scores ward 0.987 at the same compression ratio, because it routes hierarchically and its leaves stay pure. Raising threshold does not help — flat from 0.5 to 14 and worse above, since near-orthogonal sparse rows give the gate nothing to merge. Recorded rather than papered over; the sparse path documents it and points at projection="svd".

• fit_predict_sparse returned ARI ≈ 0 wherever max_leaves was smaller than the row count, and the dense path scored 1.000 on the same matrix. Labels change on every input that spends the leader budget. Two faults, one geometry. (a) Past the budget the leader pass dropped the proximity gate entirely and absorbed each further row into its nearest leader — but for near-orthogonal sparse rows a leader of n members has ‖μ‖² ≈ ‖x‖²/n, so its distance falls toward ‖x‖² as it grows while a singleton's stays near 2‖x‖². The first leader to take a second member was thereafter nearer to every remaining row, and it swallowed the stream: 4001 of 6000 rows in one leader, 1999 singletons beside it on the examples/10_sparse_highdim.py fixture at max_leaves = 2048. A leader may now hold at most 32 × (rows / max_leaves), the same constraint the dense tree offers as balance= — soft, so a row still lands somewhere when every leader is full. The multiple was swept 1 → 256 against the uncapped path on two workloads that pull opposite ways; the table and the reason a tight cap is worse than none are in sparse.rs. (b) Each row was then labelled by its nearest micro-cluster, an argmin of ‖μ_i‖² − 2⟨x, μ_i⟩ that over L centroids of one to six sparse rows is decided by ‖μ_i‖² — by how many terms those rows happened to carry — rather than by overlap. The centre-based heads (kmeans / xmeans / spherical-kmeans) now label each row by its nearest cluster centroid, which is the partition the head defines and the rule the dense estimator has used since 0.6.0; it is also O(nnz·k) against the old O(nnz·L). Neither fix works alone — A/B'd on the shipped settings, that fixture stays at −0.000 with only (b) and reaches 0.078 with only (a); with both it is 0.976 against the dense tree's 1.000. On 20-newsgroups TF-IDF the raw sparse row goes 0.004 → 0.038 and projection="svd" at max_leaves=2048 goes 0.152 → 0.195, the best that table has carried. The heads with no centre rule are dealt with in the entry above. Tables in bench/RESULTS.md.
• Mapper dropped a microcluster outright when lo + (hi − lo) rounded below hi: the cover recomputed the last bin's upper edge from the step rather than from the observed range, so the leaf holding the maximum lens value could fall outside every bin. It happened on one seed in three of the fixture that found it. The outer bin edges are now the observed range.
• feature="fd" was reporting up to a third less scatter than the leaf holds. The Frequent-Directions shrink subtracts the lower-median squared singular value from every direction, and that mass was discarded outright — so ssd(), which the radius, the absorption gate, D2/D3 and the k-means++ potential all read, under-reported the leaf. Measured share of the true scatter still reported: 0.670 at d=64, ℓ=16, 0.660 at d=784, ℓ=16, 0.897 at d=784, ℓ=64. The discarded trace is now banked, so ssd() is exact — asserted against the Full feature in a unit test and a property test, both of which fail on revert.

Handing that trace back to the covariance shape is a modelling choice, since the shrink destroys the directions along with the magnitudes, and both candidates were measured on the gmm head. Spreading it isotropically over all d is refuted: it loses in all six cells, by up to 0.32 ARI, because at d ≫ ℓ it fills hundreds of directions the sketch never saw. The shipped completion scales the retained directions up to carry the trace, which gains +0.05…+0.09 ARI on 64-dimensional digits and costs −0.04…−0.09 on 784-dimensional MNIST — the opposite error, and the smaller one. ward is byte-identical across all variants and kmeans moves only within seed spread, so the gmm head is the only one that can see this at all. Table in bench/RESULTS.md. - A rebuild now refolds every node CF from its children instead of trusting the push chain. A node CF is defined as the merge of its subtree, but insertion maintains it by pushing each point along its leaf→root path — so the root holds one sequential chain of N merges, which is exactly the O(N·u) accumulation Higham's analysis warns about. Measured on an f32 tree of 200 000 offset points, root scatter against the f64 fold of the same leaves, read right after a rebuild: 5.1e-5 before, 6.7e-8 after — a factor of 750. This is not bookkeeping: descend routes on node CFs, so their drift decides which leaf a point lands in. Labels can therefore move on float32 trees; on float64 the correction is at the last bit.

The task this closes (#44) proposed pairwise accumulation, and that half is refuted: Higham's O(m·u) vs O(log m · u) separation needs a long chain, and at the leaf counts this library runs at (180–4000) a sequential fold of the leaves and a balanced one measure 2.2e-7 against 3.2e-7 — indistinguishable. What was worth 750× is folding from the children at all. The halving is kept because it is free and the right shape, not because it was measured to matter.

Between rebuilds the root drifts again (5.1e-5 by the end of the same stream), so summary() documents the trade and points at folding leaf_features() when the aggregate is the answer.

• absorb was validated in two places that had to agree and did not. The Betula estimator restated the accepted names and the chi2_scale > 0 rule inline, next to resolve_gate's own copy of both. Adding a gate to one left the other rejecting it — precisely the drift the shared ABSORB_CHOICES constant exists to prevent, which it could not, because it shared only the message and not the set. The eager check now delegates to resolve_gate and discards the gate it builds, so the names and the rule have exactly one owner. No behaviour change for existing values.
• method="hdbscan" counted min_samples and min_cluster_size in leaves, not in points. hdbscan.rs thresholded a leaf count while its stability term used point mass, so the two arguments meant one thing on a summary and another on one-feature-per-point, and nothing said so. The failure that surfaced it: at N = 500 000, max_leaves = 2000 and the point-level min_cluster_size = 1250 a sklearn.cluster.HDBSCAN user would transfer unchanged, the head asked for 1 250 of 2 000 leaves, found no admissible split anywhere, and returned every point as noise with no warning — ARI 0.000. The threshold was also not scale-free, changing meaning whenever max_leaves changed, which is a tuning knob rather than a property of the data.

Both arguments now count points: the cluster-size gate reads the same node_mass the stability term already used, and the core distance is the smallest radius enclosing min_samples points' worth of weight rather than the min_samples-th nearest feature. On unit weights — one feature per point — every quantity is an integer count and the behaviour is byte-identical, so nothing changes for threshold=0.0, max_leaves >= N. On a summary the labels change, and for the better: the n = 500 000 blobs benchmark row moves 0.000 → 0.478, which is exactly its n = 100 000 value — the same question now gets the same answer at both scales.

A second finding, measured but not fixed here (task #72): min_samples below a single leaf's mass is enclosed at radius zero, which collapses every core distance and degenerates HDBSCAN* to single linkage. Set it above N / max_leaves; docs/USAGE.md now carries the curve and bench/RESULTS.md the comparison against fast_hdbscan that exposed it. - distance="average" panicked with an index-out-of-bounds instead of clustering. A node split picks the farthest pair of children as seeds and then assigns every child to the nearer seed — which silently assumes between(cf, cf) == 0. That holds for euclidean, manhattan and ward, but not for average (D2): the average inter-point distance between a cluster and itself is 2·S/n, not zero. A seed with enough of its own scatter therefore measured further from itself than from the other seed and joined the wrong group; once every child did that, the sibling node was created with no children, and the next descend indexed children[0] on it. Reproduced on load_digits at max_leaves ≥ 800 with any feature.

The seeds now anchor their own groups by construction, so both sides are non-empty whatever the routing measure does — the same root cause the seed scan already guards against by skipping i == j, one step later in the same function. This is the second half of that fix, and it is what the non-metric CF distances of the planned linkage and Bregman work need in order to be safe. Labels are unchanged on every metric route (48/48 (dataset × feature × head × distance × seed) label arrays hash identically across the fix), so nothing published moves.

Fixed

• The head table said method="kmeans" needs a k. It does not, and has not for some time. n_clusters=0 runs an X-means BIC selection over the leaves, exactly as the table's own prose two lines below it already claimed for "the parametric heads". Nothing caught the contradiction because the wrapper's auto-k test enumerated gmm / gmm-full / mppca / ward and omitted kmeans; it now covers kmeans too, and fails on a build with the auto arm removed. Behaviour is unchanged — only the documentation was wrong. The xmeans doc comment claimed its BIC is minimised while the code maximises it, which is the opposite of the convention the five other BIC sites in the crate use; it now says so.

• The BIC behind method="kmeans" with n_clusters=0 dropped the within-leaf scatter, and would answer k = n_leaves on any input allowed to look that far. Labels change. The score summed only the between-leaf part, Σ_l n_l‖μ_l − c‖², under the reasoning that the within-leaf scatter Σ_l S_l is the same for every k and therefore cancels. It is constant, but it does not cancel: it sits inside ln σ̂², so the two scores differ by −½·nr·d·ln((W_k − C)/W_k), which depends on k through W_k. Worse, as k reaches the leaf count the between-leaf part goes to zero, σ̂² hits its 1e-12 floor and the BIC diverges to +∞ — one cluster per leaf beats every real answer. The correct quantity is the full Σ_l (S_l + n_l‖μ_l − c‖²), the CF identity, which is exact for any whole-leaf partition and is already what KMeans::inertia accumulates. Nothing caught it because AUTO_K_MAX = 20 capped every shipped path far below the leaf count and the three tests covering the sweep all ran at k_max ≤ 8 against 100+ leaves; it surfaced only once a head whose cap is the leaf count went looking. On a 24-leaf, 6-blob fixture in 10-D the sweep selected 24 before and 6 after, and the regression test fails on a build with the term removed.

Notes

• CURE α-shrinkage (Guha, Rastogi & Shim 1998) was implemented, measured and refuted; not shipped. The full head — farthest-first representatives per cluster, shrunk toward the cluster centroid by α, linked by the closest pair — dropped straight onto the Anderberg driver, and its anchors hold (α = 1 reproduces UPGMC exactly at every k). It still does not earn a place. Against the heads already here at max_leaves=600, median of seeds 0/1/2: moons (noise 0.06) 1.000 in 0.018 s, which hdbscan also reads in 0.013 s without needing k or α; circles 1.000, same story; and then moons at noise 0.10 — the same generator, one step noisier — 0.039 against ward's 0.516 and spectral's 0.733, and 5 % uniform noise over blobs 0.008 against ward's 0.691. It never beats a shipped head on any fixture measured.

α is the reason. On moons-.06 it reads 0.032 at α = 0.05, 1.000 at α = 0.10 and 0.044 at α = 0.15; the window moves with the representative count (α = 0.05: 0.032 at c = 10, 1.000 at c = 20, 0.000 at c = 40) and with the leaf budget (0.983 at 300 leaves, 0.032 at 600, 0.964 at 1200). A setting that narrow and that unstable can only be found by scoring against the labels being computed.

The mechanism is a fact about cluster features rather than about CURE, and it is the part worth keeping. The dendrogram is right and the cut is wrong: at α = 0.05 and k = 2 the head returns clusters of 5999 and 1 points, while k = 3 returns 3000 / 2999 / 1 at ARI 0.982. That is single linkage's outlier failure, and shrinkage is precisely Guha's cure for it — but on a summary an outlier is a low-mass leaf, so the linkage meets it as a cluster of one member with one representative sitting on its own centroid, where (1 − α)·μ + α·μ = μ for every α. Shrinkage damps an outlier inside a cluster and can never damp an outlier that is a cluster. The paper's own outlier-elimination stage, transferred to leaf mass, rescues the clean fixtures and leaves the noisy ones ranked backwards: on blobs + noise the best cell is α = 1.0 (0.794), which is plain centroid linkage with no CURE in it at all. Any future rim-based head here needs a mass-based answer to that, not a geometric one.

[0.7.0] — 2026-08-23

Never published, and the heading below is deliberately not a link. The version was bumped and this section written, but no v0.7.0 tag was ever created, so nothing here reached PyPI, crates.io or Zenodo and there is no commit range to compare against. Everything in it ships in 0.8.0, whose compare link spans v0.6.0...v0.8.0 and therefore covers both. The section is kept as the record of what that heading described.

Changed

• Rust 2024 edition; minimum supported Rust version 1.82 → 1.85. The migration needed two source changes in total (a binding mode made explicit in gmm_toeplitz.rs); the rest is rustfmt's 2024 style edition collapsing short if/else onto one line. 1.85 is the edition's own floor.
• A rebuild no longer rebalances the whole tree every time it compacts. 0.6.0 set the leaf count exactly by merging the k closest sibling pairs, which is also the policy that grows the absorption threshold by the smallest amount that reaches the count — so the 10 % headroom refilled almost at once and the rebuild fired again. Each one then called the rebalance pass, which clears the node and entry arenas and re-routes every surviving entry: on 1 M blobs at max_leaves=4000 that is 385 rebuilds × ~3600 entries ≈ 1.4 M descents on top of the 1 M point inserts, and it made fit ~3.2× slower than 0.5.0 on the tree-dominated path. A profile hides it, because the rebuild's cost is descend and try_absorb.

Compaction merges strictly inside a leaf, so mass is conserved per node and the tree is already valid when it finishes — the rebalance is a re-partitioning of leaves that mix two clusters, not a repair, and the drift it corrects accumulates with merges rather than with compactions. It now runs once the merges since the last rebalance reach max_leaves. Measured on 1 M blobs at max_leaves=4000, median of 5 fits, arms alternated on an idle machine: 0.743 s → 0.230 s, 385 rebuilds → 31, leaf utilization 92.4 % → 92.6 %, ARI 0.8597 → 0.8598. Labels move wherever a rebuild actually fires — three seeds each, median [min–max]:

	0.6.0	this release
mnist-20k kmeans	0.2275 [0.1719–0.2510]	0.2571 [0.2275–0.3263]
mnist-20k gmm	0.2842 [0.2547–0.2842]	0.2788 [0.2668–0.2788]
mnist-20k ward	0.3272	0.3666
covtype-20k kmeans	0.0528 [0.0506–0.1061]	0.0878 [0.0755–0.0999]
covtype-20k gmm	0.0843 [0.0782–0.0962]	0.0882 [0.0820–0.1129]
covtype-20k ward	0.0910	0.0910

digits does not move: at n = 1797 against a 4000-leaf budget it never rebuilds. Widening the compaction margin to max_leaves / 6 was measured as the alternative and rejected — it reaches only 0.369 s and takes covtype utilization down to 83.7 %, which is the 90–99 % band the exact-count compaction was written to provide. CFTree gains one usize field, #[serde(default)], so persistence snapshots written by 0.6.0 still load and SCHEMA_VERSION stays at 2 — verified by a snapshot from the released 0.6.0 wheel, committed at tests/data/v2_0.6.0.betula alongside the attributes it should reproduce. test_save_load_roundtrip writes and reads with the same build and cannot see a schema break at all; dropping the serde(default) makes the new test fail with missing field merged_since_rebalance, which is the claim itself.

Fixed

• The vmf head reported the same concentration for every tight cluster in high dimension. κ was clamped at 10⁴, which in 784 dimensions binds from R̄ ≈ 0.9616 — ordinary for embedding data. At R̄ = 0.99 the clamp understated the Banerjee estimate 3.9×, and at R̄ = 0.999 39× (true κ ≈ 3.91·10⁵, reported 10⁴). Every saturated component therefore carried an identical, wrong κ, which biases movmf_auto's BIC uniformly and can move the k it selects.

The comment blamed numerical stability; the cause was cost. log I_ν(κ) came from an ascending power series whose peak term sits at m ≈ κ/2, so it is O(κ) inside the EM loop, and its own m > 200_000 stop truncates before the peak above κ ≈ 4·10⁵ — raising the cap alone would have traded a wrong answer for a slow one and then for a wrong one again. log_iv now branches at κ = 10⁴ to DLMF 10.41.3, the uniform asymptotic expansion for large order, which is O(1) and correctly rounded there — measured against 50-digit arithmetic over ν ∈ [1, 2047], κ ∈ [10⁴, 10⁶], the f64 result is within 0.8 ulp. The cap rises to 10⁶, set by representability rather than by the normalizer. The expansion divides by ν = d/2 − 1 and so does not exist below d = 4; there the series remains the only evaluator and the cap stays at 10⁴.

Labels move for vmf in high dimension, and n_clusters=0 can select a different component count. No published benchmark table uses the head, so no committed number changes. - predict could not return a cluster whose radius is zero. The Voronoi rule that predict applies for the centroid heads (kmeans, spherical-kmeans, ward, …) filters out clusters with no weight, but the stats helper it reads yields radii before weights for clusters — the opposite of the order it uses for leaves — and the rule destructured it in leaf order. Every degenerate cluster (a singleton, or one whose members coincide) was therefore dropped from the rule and became unreachable: the point sitting exactly on such a centre was labelled with some other cluster, at a distance three orders of magnitude larger. Only the filter was wrong; the centres themselves, cluster_sizes_ and cluster_radii_ were always correct. - Leiden was not reproducible: the same graph and the same seed gave different communities on successive calls in one process. refine (and, on a tied graph, one_level) picked its target sub-community by scanning a HashMap, and kept whichever tied candidate the iterator produced first. std::collections::HashMap draws a fresh RandomState seed per instance, so the candidate order — and with it the tie-break — changed between calls; seed only ever controlled the node visit order. aggregate had the same problem one level down: it emitted each super-node's adjacency row in HashMap order, and since floating-point addition is not associative, the weighted degree summed from that row differed run to run, moving every gain computed from it. Both now scan in community / neighbour index order. Labels change for graphs that carry tied gains — symmetric ones above all — and are now stable. Measured: reverting the three hunks in place and re-running every betula-leiden row of bench/results_quality.csv and bench/results_real.csv against the same build reproduces all nine to the last digit, so no published benchmark number moves with this fix. - Documentation shipped with 0.6.0 still described 0.5.0. The verified-suite line in README.md and DESIGN.md claimed 185 Rust + a 213-case Python suite; the tree now carries 404 Rust unit + 4 integration tests (402 + 4 with --no-default-features, 8 more for the cli binary) and a 230-case Python suite (229 passed / 1 skipped without the optional scipy / networkx test dependencies). DESIGN.md also still described predict_proba as the per-leaf responsibility matrix, which 0.6.0 replaced with scoring the row itself under the fitted mixture, and its crate layout predated src/mixture.rs. docs/api.md omitted consensus / ConsensusResult although both are exported from __all__; SECURITY.md still named 0.1.x as the supported line. - Every published benchmark number is re-measured, and several move down. bench/RESULTS.md, all eight results_*.csv, the four plots and the README headline were a single run of 2026-07-18 against a 0.2.0 build. They are now the median of seeds 0, 1, 2 against this tree, with per-cell min/median/max in new results_*_spread.csv sidecars and the seed list in results_seeds.json — because on the synthetic sets every row moves by more than 0.05 ARI across three seeds, and the old single-run page had in several places published the top of that range. Corrections worth naming: digits k-means is a tie (0.467 vs scikit-learn's 0.468, not 0.568 vs 0.468); the spectral head matches SpectralClustering at 1.0–1.5× its speed, not 3–5× (betula's time barely moved, scikit-learn's fell threefold between versions); normalize=True on MNIST is now a wash — its 0.203 baseline is gone and the off-vs-on sign flips between seeds — though it still earns its place on digits (k-means 0.467 → 0.569); and CF-weighted NMF is slower than scikit-learn's at every N measured, its win being the ±0.00 seed spread against ±0.37. The covtype loss to sklearn-birch is confirmed as a loss on the merits: matching Birch's compression ratio leaves its ARI unchanged (0.132 vs 0.131) and handing betula Birch's 11 774 leaves makes every head worse. The 23 notebooks still carry 0.2.0-era outputs and their provenance banners still say so. - References to ../../math_improove/ in DESIGN.md, src/distance.rs, src/feature.rs and src/stats.rs pointed outside the repository; they now say the working notes are local-only. - CI on macOS and Windows broke the moment .python-version was added. That job builds a deliberate 3.12 environment (uv venv --python 3.12) and installs its tools into it; with a repo-wide 3.14 pin present, a bare uv run no longer accepts that environment, recreates it, and discards the tools — error: Failed to spawn: maturin. Both uv run calls there, and the two in the mutmut job that would have hit the same wall on its next schedule, now pass --no-sync. - Mutation testing had never once produced a result, and covered a third of the crate. All seven scheduled runs since 2026-07-06 were killed by the 90-minute cap: measured on the last of them, 90 minutes at -j2 completed 98 of 1011 mutants (9.7%). Sharding alone does not fix it — a 12-way split still leaves ~115 mutants per job against a suite where a single method-dispatch test costs 25 of its 27 seconds, and every mutant pays it in full. --test-tool nextest does, because it stops at the first failing test: measured A/B over the 13 mutants of src/kernels.rs at -j4, cargo test took 525 s and reported 8 caught with 2 timeouts, nextest took 144 s and reported 10 caught with none — 3.6× faster and more accurate, since both former timeouts resolve to caught. That makes full coverage affordable, so examine_globs is gone: the scope is now the whole crate, 4310 mutants (still excluding python.rs and bin/, which have no Rust-level tests and would only yield false survivors), split 24 ways at ~180 mutants per job. A summary job collates every shard and fails loudly on a survivor. .cargo/mutants.toml records the per-module mutant counts and the measured throughput so the matrix can be resized from numbers rather than guesses. - Four mutants had been surviving in src/distance.rs unread inside a cancelled run's artifact. Radius::between was pinned only by a case where na·nb and na/nb coincide (2·1) and where one scatter was zero, so *→/ and +→- both went unnoticed; MahalanobisChi2::maha_sq was only ever exercised at mean zero, where x − μ and x + μ agree, so the sign of the residual was unconstrained. Two tests now pin both. The fourth (||→&& in VarianceIncrease::between) was an equivalent mutant — with exactly one weight zero the formula already returns zero — so the guard is now written on the sum, nab = na + nb, matching Radius::between directly below it. Same result for every input; the operator is simply no longer free. - The mutation job failed whenever any mutant survived, which meant it never carried a signal. That sounds like the strict setting, and it is the useless one: the crate has never had an empty survivor set — some survivors are provably equivalent mutants and no test can kill those — so the job was red every single week and a new hole was indistinguishable from the standing debt, which is exactly how the four src/distance.rs survivors above went unnoticed from July to August. It now fails on survivors that are not in a committed mutants-baseline.txt. That file records the set keyed file:line:col, grouped by module, each group carrying either the argument for why no test can close it or the constraint a killing fixture would have to satisfy — a claim of equivalence without an argument does not belong in it. The summary job reports both comm directions, so an entry that has since been killed surfaces as "trim this" rather than quietly padding the debt. Two failure modes of the job itself are closed with it: a shard that uploads no artifact is now an error rather than a smaller total that reads as good news (run 32347677860 lost shard 39 after 88 minutes against a ~20-minute median, and 95 of 96 collated into a number that looked complete), and the matrix is capped at max-parallel: 12, below the account's 20-slot limit, after the uncapped 96-shard run left ordinary CI of the same commit queued for two hours. - Reference and boundary fixtures close 48 of the recorded survivors — mutants-baseline.txt goes from 361 entries to 313, with src/tree.rs at 22 → 2, src/clustering/vmf.rs at 38 → 19, src/clustering/gmm.rs at 19 → 13, src/clustering/nmf.rs at 15 → 13 and src/clustering/kprototypes.rs at 4 → 3. Each fixture is verified the only way that means anything: the mutation is applied at its own line:col, with an assertion that the line still holds what the entry says it does, and the test is confirmed to fail. Four of them closed mechanisms nothing had constrained. randomized_svd had no reference check at all — the range finder is deliberately insensitive to its own sketch, so an end-to-end NMF assertion cannot see whether the sketch or the power iterations are computed correctly — and is now compared against a dense symmetric eigendecomposition of XᵀX and the Eckart–Young optimum it yields. A rebuild's grown absorption threshold, widest · (1 + 4ε), decides absorption for the whole remainder of a run and was pinned by nothing. And split seeded both its sides from a scan that, one +→* away, compares a child with itself: harmless under CentroidEuclidean, where a self-distance is zero, but Radius and AverageIntercluster both carry a child's own scatter into between, so a high-scatter child is further from itself than any two children are from each other and the split degenerates into a size tie-break that separates coincident entries. spherical_lloyd had the same shape of hole from the other side: its independent reference carries the same if !changed && it > 0 { break; } line, so no comparison against it can see that guard at all — what can is the Lloyd fixed-point condition itself, asserted on a fixture slow enough to still be moving at the second assignment. The Rust suite is 399 lib + 4 integration tests as a result.

0.6.0 — 2026-08-19

Fixed

• predict returned the label of an approximately-found leaf, not the label of the model's own partition. k-means assigns a point to its nearest centre; the labels came instead from routing the point down the CF-tree to a leaf and reading that leaf's cluster. The descent is greedy, so it is an approximate nearest-microcluster search, and in high dimension it is wrong often: measured disagreement with the model's partition of 2.7% (digits), 14.9% (MNIST-20k) and 18.1% (20-newsgroups TF-IDF) of points, and on TF-IDF only 14 of the 20 clusters the head found were reachable by descent at all — the estimator reported n_clusters_ = 20 while predict could never emit six of them. The argmin over k centres is also the cheaper of the two at k ≪ n_leaves. Measured ARI, max_leaves=2000, median of 3 seeds: MNIST-20k kmeans 0.397 → 0.417, 20-newsgroups kmeans 0.050 → 0.060 and spherical-kmeans 0.071 → 0.099, digits kmeans 0.667 → 0.663. Applies to method="kmeans" and "spherical-kmeans", whose objective is "assign to the nearest centre" (the spherical head compares unit-normalized centres, where the Euclidean argmin and the cosine argmax agree). Every other head keeps the microcluster route: the mixture heads (GMM, Toeplitz-GMM, movMF) assign by maximum posterior — which weighs each component by its own covariance and mixing weight, a different partition rather than a slower route to the same one — and Ward / Spectral / Leiden / HDBSCAN / scale-space clusters need not be convex at all. Constrained runs (fit_constrained) also keep it, since a Voronoi rule is free to violate the pairwise constraints COP-KMeans satisfied. Measured separately: a nearest-centre rule would raise the GMM heads too (digits 0.504 → 0.576, MNIST-20k 0.343 → 0.418, 20news 0.053 → 0.073), which makes a true max-posterior predict for them the obvious follow-up — see the next entry, which ships it.
• The mixture heads now label a point by maximum posterior instead of by tree descent. gmm, gmm-full, vmf and gmm-toeplitz{,-full,-gs} score the point itself under the mixture they converged to — argmax_c [ln π_c + ln p(x | θ_c)] — rather than routing it to a leaf and copying that leaf's label. Each head hands its converged parameters to a Mixture (diagonal / full-Cholesky / stationary / vMF kernels), so the density that labels a point is the density that was fitted: the floored diagonal variances, the ridge-regularized Cholesky, the AR predictor bank. predict_proba normalizes the same scores, so predict_proba(X).argmax(1) == predict(X) by construction, and it now scores the row rather than the row's microcluster. A component that no leaf claims is silenced, so predict can only emit a label the fitted partition uses. Measured ARI, median of 3 seeds, defaults (max_leaves=2000), old rule → new rule: 20-newsgroups TF-IDF gmm 0.027 → 0.054 (72.1% relabelled) and vmf 0.021 → 0.026, digits gmm 0.489 → 0.507, gmm-full 0.722 → 0.754, vmf 0.643 → 0.631, blobs gmm 1.000 → 1.000. It costs accuracy on raw image pixels at a fine leaf budget, and that is not hidden here: MNIST-20k gmm 0.340 → 0.185 (52.1% relabelled). The cause is the model, not the rule: on the same fit a nearest-centre rule, which ignores the covariance entirely, scores 0.378 — above both — so it is the fitted covariance that costs the ARI. A diagonal covariance treats 784 correlated pixels as independent, so scoring a raw point sums 784 separate penalties and the per-dimension modelling error accumulates; the leaf-level score damps this through the −½ tr(Σ_c⁻¹ Σ_i) term, which a single observation genuinely does not have, and the coarser the leaves the more that term smooths. The loss is accordingly budget-specific rather than head-specific: at max_leaves=300 the new rule wins for gmm (0.206 → 0.239) and for gmm-full (0.207 → 0.212) alike. For raw images prefer kmeans or a projection.
• The diagonal GMM could drive 1/σ² to 1e12 in a dimension with no spread. The variance floor was expressed purely relative to the spread of the dimension itself (1e-3 · gvar_d + 1e-12), so a constant or near-constant dimension — an always-zero border pixel, a term absent from every document — got a floor of essentially zero and its precision was bounded only by 1e-12. Compounding it, global_variance summed only the between-leaf term while its own doc comment claimed between + within, so the scale feeding both the floor and the NIG prior was under-estimated in proportion to how hard the tree had compressed. Neither showed while the density was only ever evaluated at leaf means, which sit close to any component mean; a raw observation does not. The floor now also carries a global term (VAR_FLOOR_REL · mean_d gvar_d, mirroring the full head's ridge = 1e-6 · tr(gcov)/d) and global_variance adds the per-leaf term, as global_cov always did. Measured ARI, median of 3 seeds, before → after, as leaf-descent / max-posterior: MNIST-20k gmm 0.317 → 0.340 / 0.135 → 0.176, 20-newsgroups gmm 0.029 → 0.027 / 0.051 → 0.059, digits gmm 0.471 → 0.479 / 0.484 → 0.496. gmm-full, vmf and blobs are byte-identical, as they must be — they never used this floor. The fix is real but partial: it lifts the MNIST point-level ARI by 30% relative and still leaves it below leaf descent, for the modelling reason given in the entry above.
• Both GMM heads seeded every component with the global covariance, so a junk dimension could decide the first E-step. gmm started each component's variance at max(gvar_d, floor_d) and gmm-full at gcov — the spread of the whole dataset. In a dimension that separates the clusters that spread is dominated by the between-cluster distance and says little about the within-component scatter: on four blobs eight units apart it overstates it 34×. A dimension whose spread is already purely within-component — a sparse binary feature, a near-constant pixel — gets no such inflation, so a uniform-looking seed is in fact a skewed metric that weighs the junk dimension up to 889× more per dimension, and one spike there costs more (51.7) than the entire blob separation (3.9). Both heads then converge with every component mean sitting on the global mean. Reproduced with four trivially separable blobs plus twelve binary columns that are 1 with probability 0.02: ARI 0.000 → 1.000 (gmm) and 0.001 → 0.985 (gmm-full), while kmeans, ward and hdbscan scored 1.000 on the very same CF-tree throughout — as did scikit-learn's GaussianMixture(covariance_type="diag"), which seeds from its k-means partition. Each component is now seeded from its own k-means cluster (the same within-cluster accumulation the M-step performs, the same floor and ridge, the global value as the fallback for an empty cluster). On real data this lands EM in a different local optimum and the effect is small in both directions (max-posterior ARI, median of 3 seeds): digits gmm 0.496 → 0.507 and MNIST-20k gmm 0.176 → 0.185, against digits gmm-full 0.777 → 0.754 and 20-newsgroups gmm 0.059 → 0.054. vmf never had the defect — it seeds κ from the hard assignment — and gmm-toeplitz starts from random responsibilities; both are unchanged to the last digit, as are all non-mixture heads.
• CF-tree: the rebuild spent a third of the leaf budget, and collapsed the tree outright in high dimension. max_leaves is a resolution budget — the summary handed to the global clustering is only as fine as the leaves actually kept — but the rebuild grew the threshold to the mean within-leaf nearest-sibling gap and then merged whatever fell under it, which is a prediction, not a target. Measured utilization of the budget: 65.6–96.2% (mean 80.6%) across digits, blobs at n = 20 000/100 000 and d = 10/50, and uniform d = 20. On 3000-dimensional TF-IDF it failed catastrophically instead — 3 leaves against a 2000 budget, one of them holding 85.7% of the mass, ARI 0.0 for every method and every seed. The cause is concentration of measure: the achievable leaf count is near-discontinuous in the threshold (7755 leaves at threshold=1.0, 12 at 1.3), so no threshold-first policy has a safe value there. The rebuild now merges the k closest sibling pairs with k set by the budget and reads the grown threshold off the widest gap it took — k is exact, and merging is capped at one pair per entry, so the cliff cannot be stepped over. Utilization is now 90.0–99.0% (mean 93.5%) on the same datasets and 90.0/94.1% on TF-IDF (was 0.6/0.2%), where ARI goes 0.0 → 0.071 with the spherical-kmeans head.
• CF-tree: the rebuild conflated reducing the entry count with rebalancing the node structure. Merging two entries inside their own leaf node leaves every node CF in the tree exactly unchanged — a node's CF is the merge of its subtree, and merging two of its children does not change that multiset union — so mass is conserved per node, no ancestor needs touching, and no leaf can be emptied. The count reduction is therefore done in place, in one O(Σ_leaf child_count²) scan plus an O(m log m) sort, and the reinsertion pass that follows now merges nothing: it only re-routes entries so that leaves re-partition around the geometry the data actually has, which compaction alone cannot do (it can shrink a leaf that mixes two clusters but never split it). Keeping absorption on during that pass is what walks off the concentration cliff — measured, it collapses a d = 50 blob mixture to 9 leaves against a 500 budget. Rebuilds became more frequent (each shaves 10% rather than overshooting) and individually cheaper; end-to-end build time moved −6% to +98% across the eight configurations, with the two largest (blobs n = 100 000, d = 50) essentially unchanged at +1% and +3%. n_rebuilds_ counts a different, cheaper unit of work than before and is not comparable across the change.

Quality follows the head, not the leaf count, and the finer summary is not uniformly better: on digits (5 seeds, median) the kmeans head goes 0.646 → 0.667 while the gmm head goes 0.593 → 0.458 at max_leaves=500 — the diagonal-covariance head degrades as leaves get smaller and is best on that dataset at ~230 leaves, which is now a max_leaves choice rather than something the rebuild does behind the caller's back. Blob mixtures stay at ARI 1.000 throughout. - CF-weighted NMF: the solver stopped after 5 sweeps regardless of max_iter. The convergence test compared |prev − err| against tol · prev with prev seeded at +inf; IEEE-754 says inf <= inf, so the first check always fired and the iteration budget was dead code. Every other EM loop in the crate guards this with it > 0; nmf.rs was the only one that did not. The budget is now honoured, and a regression test asserts the residual keeps falling as the budget grows. - CF-weighted NMF: most components collapsed to zero and never recovered. Two independent causes, both fatal because zero is an absorbing state for HALS and for the multiplicative updates. (1) The randomized SVD recovers the right vector as v = Bᵀu/σ; below the noise floor that division amplified round-off into a vector of arbitrary magnitude, so a rank-deficient triplet seeded a component wildly out of scale. Such triplets are now cut at the LAPACK numerical-rank threshold and reported as exact zeros. (2) NNDSVD's zero-fill used mean(X), but a filled component is a rank-1 block of constant magnitude and r of them add up, so the fill swamped the data as the rank grew. The init now uses the ar fill, mean(X)·U(0,1)/100 — small enough to stay a perturbation, random enough to break the degeneracy a constant fill creates. Measured on a rank-12 matrix at rank 32: initial relative residual 13.5 → 1.2, dead components 28/32 → 3/32 (3 is the honest count — the data is rank 12), and the converged residual 270× lower. On digits at rank 24: 15/24 → 0 dead, reconstruction error 0.33 → 0.20, downstream ARI 0.54 → 0.61. The reconstruction error now matches scikit-learn's own NMF to three decimals at every rank tested (10, 16, 24, 32), against being 2.3× worse before. - CF-weighted NMF: an arbitrary per-component scale reweighted the downstream clustering. NMF is invariant to (W D, D⁻¹H), so the optimizer left whatever split it landed on — measured spreads of 70× between component scales on a converged fit. Because the codes are consumed as a Euclidean feature vector by the Phase-3 head, that scale acts as a per-dimension weight and the head clustered along whichever component drew the largest number. The factorization is now canonical (component rows unit-L2, scale absorbed into the codes, components ordered by descending energy). Measured on a 4-topic nonnegative mixture over 8 seeds: median ARI 0.63 → 1.00, seed spread 0.37 → 0.00; against the shipped 0.5.0 behaviour, median 0.92–1.00 → 1.00 with spread 0.31 → 0.00 at N = 160 000. scikit-learn's own NMF shows the same 0.37 spread on this data, unfixed.

Added

• NNDSVDar initialization for both NMF solvers (Boutsidis & Gallopoulos, 2008), built on a self-contained randomized range finder (Halko-Martinsson-Tropp) — deterministic given the seed and far better conditioned than the previous random start, which decided the basin a non-convex coordinate descent lands in.
• Betula.components_ — the NMF parts H as (projection_dim, dim), unit-L2 rows ordered by descending energy, so a row reads directly as a topic over the input features.
• Betula.reconstruction_err_ — relative reconstruction error ‖X̃ − W H‖_F / ‖X̃‖_F of the projection over the leaf centroid matrix.
• projection_max_iter (default 100) — the factorizer's own sweep budget, independent of the head's max_iter. Previously the two shared one number, so raising the clustering budget silently paid for NMF sweeps too.
• projection="weighted-nmf" / "weighted-nmf-kl" now accept sparse CSR input; only the stored values are checked for negativity (implicit zeros are already nonnegative). Measured ARI 1.00 on a sparse topic mixture, matching the dense path.
• A convergence check for the KL solver (it previously always burned the full max_iter). The Frobenius solver's check now follows the size of the update rather than the size of the objective, compared against the first sweep's — sklearn's rule. A relative test on the residual provably never fired: HALS converges sublinearly and keeps buying more than tol of relative improvement for hundreds of sweeps, so max_iter was the only brake (measured: a 4× budget cost 6× the time). The movement falls out of the sweep for free, so the check is also cheaper than the residual it replaced.

Changed

• WᵀX, the NMF sweep's hot loop, restructured for locality: measured 2.4-3.4× faster. Evaluating it per output cell (Σ_j w[j][k]·x[j][c]) walks a whole column of X for each of the r·d cells, striding d floats per step through a matrix far larger than L2. Accumulating into the small, cache-resident r×d output while reading each row of X once, sequentially, computes the same product with the access pattern the hardware wants (verified equal to 1e-9 by a unit test).
• The NMF rank is capped by both centroid-matrix dimensions (min(rank, d, M)): rank > M makes the factorization rank-deficient by construction and leaves whole components with nothing to fit.

Notes

• A-HALS extrapolation (Ang & Gillis, 2019) was implemented, measured and rejected. It lost to plain HALS at every sweep budget on two planted problems (up to 445× worse residual at 200 sweeps): the accept/reject test for the extrapolation parameter has to judge the objective at the feasible iterate, but the Gram factors on hand belong to the extrapolated point, so making it honest costs the extra O(Mdr) product the acceleration was meant to save. The cheap exact convergence check that came out of the attempt was kept.

0.5.0 — 2026-07-18

Added

• method="gmm-toeplitz-gs" — full-order Gohberg-Semencul MLE Toeplitz-precision GMM for ordered stationary signals: a Yule-Walker (Levinson) warm start at order ≤ 16, then coordinate ascent of the exact log-likelihood over the reflection coefficients (positive-definite by the |k| < 1 constraint) — the likelihood-optimal general precision (arXiv:2311.14995), completing the three-rung Toeplitz ladder of docs/adr/001-gmm-toeplitz.md. Reuses the exact GS precision E-step; BIC auto-k; a true predict_proba. Measured: competitive with the AR head on AR signals and recovers mid-lag echo structure (lags 11–16, beyond the banded cap w_max = 10) the banded head is blind to, while the dense-covariance gmm-toeplitz-full covers arbitrarily long lags (bench/toeplitz_ar_mixture.py).
• projection="weighted-nmf-kl" — KL-divergence (I-divergence) CF-weighted NMF for count data (word counts / event tallies / Poisson observations), where the Frobenius objective (Gaussian noise) is mis-specified. Lee-Seung multiplicative updates with the shared components weighted by leaf mass — the Poisson maximum-likelihood fit. The advantage is largest where counts are sparse: measured up to +0.5 ARI over Frobenius on a Poisson-count mixture at mean rate < 0.5 (0.83 vs 0.24), narrowing to a few points as counts grow and Poisson → Gaussian (examples/17_nmf_topics.ipynb).

Notes

• NMF Phase-2 warm-start / randomized range-finder were assessed and deferred (see plans/nmf-cf-weighted.md): in the CF-compressed regime the factorization already runs over the M ≪ N leaves, so the compression — not these — is the speedup; they would add API/state for a marginal gain and are not shipped.

0.4.0 — 2026-07-18

Added

• projection="weighted-nmf" (+ projection_dim) — CF-weighted nonnegative matrix factorization as a Phase-3 reducer for nonnegative data (TF-IDF / bag-of-words / event counts / spectrogram magnitudes / histograms). Rather than factorizing the raw N×d matrix (which defeats the compression), it factorizes the M ≪ N leaf centroids weighted by their mass, X̃_j = √n_j·μ_j: by König-Huygens the full-data NMF objective equals the weighted-centroid one up to the within-leaf scatter constant, so the expensive factorization runs over the microclusters — memory-bounded, O(M·d·r) — and any head (k-means / GMM / Leiden) then clusters the nonnegative codes. The solver is a dependency-free weighted HALS (coordinate descent, Gram-reuse across sweeps); the compression, not a fast NMF, is the speedup. Available on the one-shot fit_predict and the streaming Betula estimator. Signed input is rejected (no silent shifting — use vmf / spherical-kmeans or PCA / TruncatedSVD for embeddings); dense input in this release. See docs/MATH.md.

0.3.0 — 2026-07-18

Added

• method="gmm-toeplitz-full" — general (non-AR) positive-definite Toeplitz-covariance GMM for ordered wide-sense-stationary signals whose autocovariance a low-order AR cannot capture (e.g. a long-lag echo / narrowband structure beyond the AR order). Each component covariance is the dense Toeplitz matrix built from the biased (periodogram-consistent) autocovariance — positive- semidefinite by construction, made strictly positive-definite by the within-leaf variance plus a small ridge — factored by Cholesky for an exact multivariate-Gaussian E-step. O(d²) parameters, O(d³) per component; BIC auto-k at n_clusters=0; a true posterior via predict_proba. This is the general (non-AR) rung of the Toeplitz ladder recorded in docs/adr/001-gmm-toeplitz.md. On a long-lag-echo mixture (echo lag K ∈ {16, 28, 40}, all beyond the AR order) it recovers the components (ARI 0.70 → 0.97 as the window grows) where the banded gmm-toeplitz sits at chance; on AR-generated signals the two match (bench/toeplitz_ar_mixture.py).

Changed

• gmm-toeplitz: raised the internal AR-order cap w_max 6 → 10. BIC still selects the smallest sufficient order (easy signals are bit-for-bit unchanged); higher-order / MA-like signals gain headroom before the general gmm-toeplitz-full head is needed.

0.2.0 — 2026-07-18

Added

• method="spherical-kmeans" / method="vmf" — directional clustering on the unit hypersphere for L2-normalized embeddings (CLIP / face / sentence / speaker vectors), where cosine — not Euclidean — geometry is what matters. spherical-kmeans is hard cosine assignment (argmax_c μ̂·μ_c, centers re-normalized to the sphere); vmf is a soft mixture of von Mises–Fisher distributions (EM, a true posterior for predict_proba, and BIC auto-k when n_clusters=0). Both reduce each leaf to its weighted mean (n_i, μ_i), so the cluster resultant R_c = Σ n_i μ_i stays exactly mergeable — the BETULA property carries through to the sphere — and the within-leaf spread ‖μ_i‖ feeds the concentration κ (Banerjee et al. 2005), so microcluster compression does not over-estimate it. The engine L2-normalizes input automatically for these methods (get_params stays verbatim). The κ normalizer uses a dependency-free, numerically stable log I_ν(κ) (log-space series) — no Bessel library, the crate stays NumPy-only. Available on the dense one-shot / streaming estimator, the sparse (O(nnz)) path, and as the spherical_kmeans / movmf / movmf_auto Rust functions.
• covariance_weight (method="leiden" / "leiden-cpm", feature="full") — covariance-aware community detection. β > 0 adds a log-Euclidean shape term β·‖logΣ_i − logΣ_j‖²_F to the microcluster affinity graph, so two microclusters must be close in both centroid and covariance to be neighbours — useful when clusters differ by orientation / shape (covariance descriptors, motion / time-series windows, anisotropic blobs). logΣ is computed with the in-house Jacobi eigensolver (new linalg::matrix_log) — no new dependency; β = 0 (default) is the existing centroid-only affinity, bit-for-bit unchanged.
• tangent_weight / tangent_rank (method="leiden" / "leiden-cpm", feature="full") — GeoBETULA manifold-aware community detection. γ > 0 adds a Grassmann term γ · d²_Gr(U_i, U_j) (projection distance between each microcluster's rank-tangent_rank principal subspace) to the affinity, so communities must agree in centroid, covariance and manifold orientation — separating crossing / adjacent manifolds that share a centroid neighbourhood. Reuses the in-house Jacobi eigensolver (no new dependency); γ = 0 (default) leaves the graph unchanged.
• method="scale-space" — scale-space (Morse-persistence) density-mode clustering. Treats the microclusters as a weighted point set and clusters the modes of the KDE ρ_h(x) = Σ_j n_j exp(−‖x−μ_j‖²/2h²); it sweeps the bandwidth h and keeps the labelling at the most persistent mode count (the widest plateau of the modes-vs-log h curve), so it needs no k and no bandwidth and finds non-convex, arbitrary-count structure. A prominence-based mode merge (collapse peaks separated by only a shallow density valley) cleans the mode-count curve, so it is robust from 2 to ~8+ well-separated clusters and on unequal densities. Pure-Rust mean-shift over the M ≪ N leaves — cost bounded by the leaf budget, not N.
• method="gmm-toeplitz" — AR / Toeplitz-structured GMM for ordered, wide-sense-stationary signals (fixed-length time-series windows, trajectories, sensor / audio / vibration waveforms). Each component's covariance is an AR(w) process: the pooled unbiased (covariance-method) autocovariance is mapped by Levinson-Durbin to the exact Gohberg-Semencul precision Γ = (1/σ²)(BBᵀ − ZZᵀ), evaluated by the prediction-error decomposition so the w boundary positions are modelled exactly — positive- definite by construction (the reflection-coefficient clamp is the GS box constraint), O(w) parameters, order w chosen by BIC — so it stays well-posed in the N_k ≪ d regime where full covariance is singular and a diagonal model is blind to neighbour correlation. Reuses the CF scatter (no new tree machinery); a scalar stationary mean; BIC auto-k at n_clusters=0; a true posterior via predict_proba; parallel EM restarts (parallel feature). For ordered coordinates only — on generic embeddings the Toeplitz prior is wrong (use gmm / gmm-full). Based on the Gohberg-Semencul Toeplitz-precision estimator of arXiv:2311.14995; design and validation in docs/adr/001-gmm-toeplitz.md.

Fixed

• High-dimensional GMM regularization — the expected-log E-step adds a within-leaf correction (−½ Σ_d (Σ_i)_dd/σ²_kd for diagonal, −½ tr(Σ_k⁻¹ Σ_i) for full covariance) that turns over-confident when a component's own covariance goes near-singular along a low-variance direction — which is the norm in high dimensions with few effective microclusters per component. Two floors now keep the component covariances well-conditioned:
• method="gmm" (diagonal): per-dimension variance floor raised from 1e-6·gvar_d to 1e-3·gvar_d. digits (64-D) ARI 0.372 → 0.396, now ahead of scikit-learn's GaussianMixture(covariance_type="diag") (0.324).
• method="gmm-full" (full covariance): added a per-dimension floor on each component's covariance diagonal at 1e-3·gcov_dd (off-diagonals — orientation — untouched). Previously a component could be starved to zero responsibility and the recovered count dropped below k; on digits the fit collapsed to 9 clusters at ARI 0.391, and now holds all 10 at ARI 0.511 — ahead of scikit-learn's GaussianMixture(covariance_type="full") (0.402).

The floors are relative to the per-dimension global variance (not the global mean scale, which is inflated by between-cluster separation and would over-regularize tight clusters), so low-dimensional and anisotropic fits are unchanged (well-separated blobs still ARI 1.00; the rotated-anisotropic 2-D case still ties GaussianMixture at 0.887). No API change.

0.1.5 — 2026-07-04

Added

• method="leiden" / method="leiden-cpm" — graph clustering / community detection over the microcluster affinity graph via the full Leiden algorithm (Traag, Waltman & van Eck 2019): local moving → refinement (sub-communities grown from singletons along edges, so each is connected by construction — Leiden's guarantee over Louvain) → aggregation seeded from the pre-refinement partition. It discovers the community count — no k (like the density head). A resolution (γ) knob trades community count against size; the modularity objective ("leiden", γ = 1 default) has a resolution limit, the CPM objective ("leiden-cpm") is resolution-limit-free (γ on a smaller, density scale). Pure Rust — no eigensolver, NumPy-only. Best for community/blob structure at a moderate threshold; use method="spectral" for elongated manifolds. The self-tuning k-NN affinity graph is shared between the spectral and Leiden heads.
• betula_cluster.consensus(X, n_clusters, n_runs=…) — clusters X under several random insertion-order permutations and votes, turning the CF-tree's insertion-order sensitivity (Known Limitation #1) into a measurable quantity: a consensus labelling plus a per-point stability score in [0, 1] (ConsensusResult.confidence — low on unstable boundaries, high where every order agrees). NumPy-only; for the partitional heads at a fixed n_clusters. The independent runs parallelize across threads with n_jobs (the Rust core releases the GIL).
• method="spectral" — spectral clustering over the CF-tree leaf microclusters for non-convex / manifold clusters (moons, rings, spirals) that the centroid heads cannot separate. Self-tuning symmetric k-NN affinity (Zelnik-Manor & Perona local scaling), the normalized Laplacian embedding (Ng-Jordan-Weiss) via the in-house Jacobi eigensolver — no LAPACK/ARPACK, the crate stays NumPy-only — with a k-means landmark reduction above 256 microclusters so the O(m³) solve stays bounded. Dense input only; pair it with a small threshold (many leaves) so the microclusters resolve the manifold. No built-in cluster-count selection: n_clusters=0 defaults to 2.
• threshold="auto" for the Betula estimator — removes the one hyperparameter users most often have to guess. A bounded-subsample pilot fits a threshold=0 tree at the same max_leaves and reads the threshold it converges to, warm-starting the full fit near-converged instead of growing it from zero (fewer rebuild passes, lower peak leaf count on large n). Cached across refits / streaming batches; below the pilot cap it is a no-op (growing from zero is already cheap), and it is dense-only (raises on sparse input).

Changed

• Benchmarks now cover every head (spectral, Leiden added to bench/comprehensive.py) and the compression heads run at max_leaves = 4000: betula-kmeans is at exact parity with scikit-learn (blobs 0.861 = 0.861) and Ward beats raw Ward while running the full N. Docs / README / docs site surface the spectral, Leiden and consensus additions; test counts reconciled (190 Python, 158 Rust). The docs site now renders the CHANGELOG and redeploys on every published release.

0.1.4 — 2026-07-04

Added

• MapperGraph.persistence_diagram / MapperGraph.persistence(filtration=…) — 0-D persistent homology of the Mapper nerve by single-linkage union-find (elder rule, O(E log E), pure Rust). Two filtrations: "overlap" (the 1 − edge_overlap Bhattacharyya gap — a finite bar's death is the depth of a bottleneck, ranking the boolean bridges) and "lens" (the lens sublevel diagram). Essential connected-component classes carry inf death.
• Greedy weighted k-means++ init (scikit-learn's default): lower-inertia, lower-variance seeds at ~ln k× the negligible init cost over the leaves.
• objective="dbcv" for tune — Density-Based Clustering Validation (Moulavi et al. 2014, in [-1, 1]). Unlike the convex Calinski-Harabasz / Davies-Bouldin metrics (which penalise correct non-convex partitions), DBCV validates variable-density / non-convex clusters, so it is the right selection metric for the HDBSCAN-CF and DbStream density heads. NumPy-only, computed over a subsample.

Changed

• fit_predict_sparse / the _core CSR entry points now cap n_features (MAX_SPARSE_FEATURES) and validate CSR arrays through the pure-Rust sparse::validate_csr, closing an unbounded-allocation DoS where a hostile caller could force an ~8 EB allocation with a single-nonzero row.
• Docs reconciled to the current suite: 172-case Python suite, 147 Rust tests (143 unit + 4 integration under default features; the python / persistence / cli surfaces add more, 155 total).

Tests

• Mutation-testing infrastructure (cargo-mutants scoped to the CF math core, mutmut for the Python wrapper, a weekly non-blocking workflow) plus a CSR-fuzzing proptest and the two coverage gaps it surfaced (the CF-tree absorption boundary, exact tune-metric values).

0.1.3 — 2026-07-04

Added

• betula_cluster.tune — memory-aware hyperparameter search over the CF knobs, scored by an internal metric (Calinski-Harabasz / Davies-Bouldin) or ARI, with a multi-objective quality / memory / speed Pareto mode. NumPy-only by default; an optional Optuna backend (TPE / NSGA-II) via pip install 'betula-cluster[tune]'.
• Property-based tests (proptest, dev-only) for the CF-tree invariants: the clustering feature is a commutative monoid (merge is associative/commutative and equals a sequential build), folding a tree's leaf features reconstructs the whole-dataset feature, the full-covariance upper-triangular index is a bijection (incl. dim ≥ 4), and the Frequent-Directions sketch is lossless on low-rank data and never overshoots the exact scatter.
• Sparse-text benchmark (20 newsgroups, TF-IDF): the O(nnz) fit_predict_sparse path and the standard reduce-then-cluster pipeline (TruncatedSVD / NMF → k-means) vs scikit-learn, written up honestly in bench/RESULTS.md (raw high-d TF-IDF concentrates for every fast clusterer; on NMF topics betula matches sklearn).
• MapperGraph.edge_overlap — a Bhattacharyya coefficient in (0, 1] per Mapper edge, from the pooled diagonal-Gaussian summaries of the two nodes' member microclusters. Surfaced on to_networkx() edges as overlap=…, so a bridge between well-separated regions reads as a lower-weight edge than one inside a dense blob.
• Documentation site (MkDocs Material + mkdocstrings API autodoc, MathJax-rendered math) built from docs/, with a GitHub Pages deploy workflow; pip install 'betula-cluster[docs]' for the toolchain.

Changed

• Coverage floor (cargo llvm-cov, ≥95 % lines) now also measures the persistence and cli feature sets, not just the default core.
• Declared rust-version = "1.82" (MSRV) and lowered the real floor to it — the streaming heads had an implicit 1.87 dependency (u64::is_multiple_of), now rewritten. Added Documentation / Changelog project URLs.
• Docs reconciled to the current suite: 167-case Python suite, 141 Rust tests (137 unit + 4 integration), and five end-to-end use cases (README, DESIGN.md).
• Repository hardening: macOS / Windows CI test legs, an sdist install smoke test, a nightly cargo audit cron, Dependabot, and SECURITY.md / CONTRIBUTING.md / issue templates.

0.1.2 — 2026-06-28

Added

• betula_cluster.__version__, resolved from the installed package metadata.

Changed

• README repositioned: compress-then-cluster framing, the test/coverage story surfaced at the top, a "When to use it" section, and a stable-core / experimental capability split. HDBSCAN is labelled HDBSCAN-CF consistently in prose (the method="hdbscan" API string is unchanged).

Fixed

• Stale docs: the Python suite is 153 cases (was written as 123); betula-index references now point to lexindex (the indexing companion's published name).

0.1.1 — 2026-06-28

Fixed

• PyPI project description: README links to the docs, benchmarks, and examples are now absolute GitHub URLs so they resolve on the PyPI page (relative links only worked in the GitHub-rendered README).

0.1.0 — 2026-06-28

First public release.

Added

• Numerically stable BETULA clustering features (n, μ, S) (Welford/Chan updates) with four covariance models: spherical, diagonal, full (PSD via Cholesky), and a Frequent-Directions sketch (O(ℓ·d) per leaf) for very high-dimensional data.
• Memory-bounded CF-tree (Phase 1) with auto-rebuild under a max_leaves cap; optional parallel shard+merge build (n_jobs); EWMA decay for streaming concept drift.
• Global clustering heads: Hamerly-accelerated exact k-means, diagonal & full-covariance GMM-EM (expected-log E-step + NIW/MAP), Ward-HAC (nearest-neighbour chain), and HDBSCAN-on-CF; automatic cluster count at n_clusters=0 (BIC / X-means / dendrogram cut).
• χ² / Mahalanobis mass-invariant absorption gate (absorb="chi2").
• normalize=True for cosine/direction clustering of embeddings (L2-normalized rows on the unit sphere; squared-Euclidean is monotone in cosine). Doubles as the high-dimensional fix: at d≫100 raw Euclidean distances concentrate and the CF-tree collapses, but direction stays discriminative — on MNIST-784 it lifts ARI 0.04 → 0.44, beating scikit-learn (benchmarked in bench/results_real_normalize.csv). Off by default (magnitude is signal on tabular data).
• Inline auto-vectorized distance kernels (the compiler vectorizes the tight reductions per call site; target-cpu=native opts into AVX2 / AVX-512 — see .cargo/config.toml); rayon-parallel labeling.
• Python bindings: abi3 wheel (CPython 3.11+), zero-copy NumPy, float32/float64 (no upcast), GIL released during compute; one-shot fit_predict and a scikit-learn-style streaming Betula estimator (partial_fit / fit / predict / fit_predict).
• Full scikit-learn parameter protocol (get_params / set_params) — works with clone, Pipeline, and GridSearchCV. PEP 561 typed (py.typed + stubs).
• Dataset-structure inspection: microcluster_centers_/_weights_/_radii_, cluster_centers_/_radii_/_sizes_, outlier_scores, find_outliers, find_near_duplicates, near_duplicate_pairs (scored cosine pairs, exact within each leaf-block — the scalable counterpart to an O(N²) all-pairs scan), sample_representatives, assign_microclusters, summary, and n_rebuilds_ / threshold_ diagnostics.
• Mapper topological skeleton (topology::mapper → Betula.mapper() → MapperGraph): a lens (density / radius / l2norm / coordinate / eccentricity) over the microclusters, an overlapping cover, per-bin single-linkage at a data-adaptive (median-NN) scale, and a nerve graph with branch points and bridges (Tarjan); optional to_networkx(). Exploration of structure / RAG leakage / dedup, not a partition. mapper_stability() sweeps the resolution and reports the topology's persistence across scale (β₀ components, β₁ loops, branch points, bridges per resolution).
• Soft assignment & confidence: predict_proba (true posterior for the GMM heads via the per-leaf responsibility matrix microcluster_proba_; a documented centroid-distance softmax heuristic for k-means / Ward / HDBSCAN) and assignment_confidence.
• Coreset / diagnostics: export_coreset() → Coreset (leaves as weighted points — a streaming coreset), diagnostics() (compression ratio, radius percentiles, cluster mass spread), representatives(method=medoid|boundary|outlier|diverse), and cluster_profile() (JSON-able geometry for LLM cluster naming).
• memory_budget_mb: size max_leaves from a target tree-resident memory (MiB) at fit time instead of tuning it by hand; the resolved value is exposed as effective_max_leaves_.
• Drift monitoring & curation: snapshot() + Betula.compare_snapshots(before, after) (nearest-centroid match → centroid shifts / mass ratios) and active_learning_batch(strategy= "uncertain"|"outlier") (rows to review/label).
• DenStream streaming density clusterer (Cao et al., SDM 2006) over fading spherical micro-clusters built on the stable CFs (decay is centroid/radius-invariant); partial_fit / cluster / fit / fit_predict / predict (-1 = noise) + microcluster getters, sklearn-style.
• DbStream streaming DBSTREAM clusterer (Hahsler & Bolaños, 2016): fading micro-clusters connected by shared density (faded overlap mass) rather than distance, so it recovers arbitrarily-shaped clusters and keeps close-but-disconnected dense regions apart. Fixed-radius multi-assignment online; offline connects a pair when their overlap mass is ≥ alpha·min_weight. Same fading-CF core and sklearn-style API as DenStream; core::stream::DbStream in Rust.
• Streaming quantile sketches (betula-sketch, in src/sketch/): KllSketch (Karnin–Lang– Liberty, rank-error) and DdSketch (Masson et al., relative-error) — update / update_many / merge / quantile / quantiles, mergeable, bounded memory.
• Sparse input: fit / fit_predict / partial_fit / predict accept a scipy.sparse matrix (CSR-routed, rows expanded one at a time — the dense N × d matrix is never materialized). f64; this dense-tree path keeps the cancellation-free guarantee, compute O(N·d).
• O(nnz) sparse-native (fit_predict_sparse): one-shot clustering of a scipy.sparse matrix that touches only the non-zeros. Rows summarize into spherical micro-clusters keeping (n, ΣX, ‖ΣX‖², S) (so the mean, cached ‖μ‖², and centroid distance are O(nnz)) via a flat leader pass bounded by max_leaves, then a parametric head (kmeans default — robust for high-d sparse) labels each row. Uses the expanded squared-distance form, so unlike the dense path it is not cancellation-free (accurate for sparse rows far from the dense centroid); core::sparse::{summarize_sparse, nearest_sparse} is the Rust API.
• Robust insertion (huber_k): optional Huber/winsorized point updates on the streaming estimator — each point is clamped to within huber_k per-dimension standard deviations of its target microcluster before the Welford fold-in, bounding any single point's pull on the centroid (O(k·σ/n)) so stream outliers cannot stretch a centroid or inflate a radius. Off by default; zero-variance dimensions pass through and a 5-point warm-up gates the clip. The result is still a valid (n, μ, S) triple, so every downstream head is unchanged.
• Constrained clustering (must_link / cannot_link): semi-supervised COP-KMeans (Wagstaff et al., 2001) over the leaf microclusters — fit(X, must_link=..., cannot_link=...) / fit_predict(...) take (m, 2) row-index pairs. Must-link is transitively closed; cannot-link is enforced per assignment. Constraints are honoured at the microcluster granularity, so a cannot-link inside one leaf (or contradictory / over-constrained inputs) raises ValueError rather than being silently dropped. method="kmeans", dense input; core::clustering::cop_kmeans exposes the Rust API with a typed ConstraintError.
• Mixed numeric + categorical clustering (KPrototypes): k-prototypes (Huang, 1997) for mixed data. A mixed CF (MixedCf) pairs the stable numeric (n, μ, S) with a per-attribute category histogram (mode = categorical centroid); distance is ‖Δnumeric‖² + γ·(categorical mismatch), with γ defaulting to Huang's heuristic. Rows are leader-summarized into bounded mixed micro-clusters, then clustered. Standalone scikit-learn-style estimator (categorical column indices, fit/fit_predict/predict, cluster_centroids_/cluster_modes_); core::clustering::{MixedCf, kprototypes, summarize_mixed} is the Rust API.
• Command-line interface (betula, behind the cli feature): a dependency-free binary that clusters a delimited numeric file or stdin and writes one label per row; flags mirror the library (--clusters / --method / --feature / --threshold / … ; --clusters 0 auto-selects k).
• save / load + pickle (joblib-compatible) persistence (serde + CBOR via ciborium, schema-versioned).
• NaN/Inf input validation at the boundary.

Fixed

• estimate_threshold now measures the mean nearest-sibling distance within each leaf node (ELKI/BETULA-standard, O(M·capacity)) instead of a global all-pairs scan; the rebuild threshold rises monotonically (no multiplicative bump that compounded across rebuilds and collapsed the tree far below max_leaves), and rebuilds reinsert in reverse-DFS leaf order. The CF-tree build is now byte-for-byte the reference (betulars) tree shape and at speed parity with matched build flags.