diff --git a/.gitignore b/.gitignore
index d6d5142..ea1351e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,3 +10,5 @@
*.so
*.dylib
*.egg-info/
+machine.yaml
+logs/
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 5103252..2f2d6c8 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -34,7 +34,7 @@ python scripts/benchmark_suite.py \
## Release workflow
1. Update benchmark artifacts if the public performance story changed.
-2. Ensure `README.md`, the notebook, and `benchmarks/results/*.json` agree with the current implementation.
+2. Ensure `README.md` and `benchmarks/results/*.json` agree with the current implementation.
3. Tag the release as `vX.Y.Z`.
4. Push the tag to trigger `.github/workflows/release.yml`.
5. Verify the uploaded GitHub release artifacts and the published PyPI files.
diff --git a/Cargo.lock b/Cargo.lock
index 8392527..67f325d 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -165,13 +165,14 @@ checksum = "c8d4a3bb8b1e0c1050499d1815f5ab16d04f0959b233085fb31653fbfc9d98f9"
[[package]]
name = "clostera"
-version = "1.0.4"
+version = "1.0.5"
dependencies = [
"approx",
"criterion",
"ndarray",
"ndarray-linalg",
"numpy",
+ "pkg-config",
"pyo3",
"rand 0.9.4",
"rand_chacha 0.9.0",
diff --git a/Cargo.toml b/Cargo.toml
index a4eaa34..971ca15 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,10 +1,10 @@
[package]
name = "clostera"
-version = "1.0.4"
+version = "1.0.5"
edition = "2024"
[features]
-default = ["openblas-system"]
+default = []
openblas-system = ["ndarray-linalg/openblas-system"]
openblas-static = ["ndarray-linalg/openblas-static"]
python = ["dep:numpy", "dep:pyo3", "pyo3/extension-module"]
@@ -23,6 +23,9 @@ rand_chacha = "0.9"
rayon = "1.11"
thiserror = "2.0"
+[build-dependencies]
+pkg-config = "0.3"
+
[dev-dependencies]
approx = "0.5"
criterion = { version = "0.5", default-features = false, features = ["html_reports"] }
diff --git a/README.md b/README.md
index 6e768e5..5cce5c1 100644
--- a/README.md
+++ b/README.md
@@ -1,773 +1,450 @@
-# clostera: The Billion-Vector Resurrection
+
-
- 
-
+Made with ❤️ by [Synerise](https://synerise.com).
-**They told you that clustering massive high-dimensional vector collections on a single machine was a fool's errand. They said you needed a cluster, a distributed headache, and a cloud bill large enough to ruin your week. They were wrong.**
+Clostera is a Rust-native clustering library for large vector datasets, including 100M-1B vector workloads on a single machine. The public API is deliberately small: pass vectors, pass `K`, pass the metric, and either let `algorithm="auto"` choose the backend or select a concrete algorithm by name.
-`clostera` is a from-scratch Rust rebuild of the original `pqkmeans` repository, aimed at the workloads that made that project exciting in the first place: extremely large vector collections, high dimensionality, single-machine practicality, and performance that is measured rather than hoped for.
+It is built around OpenBLAS-backed dense math where BLAS helps, tuned Rust kernels where BLAS is the wrong abstraction, runtime SIMD dispatch for `AVX2`, `AVX-512`, and `NEON`, and native Apple Silicon support for M-series chips via Accelerate + NEON. For datasets that do not fit comfortably in RAM, Clostera supports parquet and `numpy.memmap` workflows so the heavy data can stay out-of-core.
-This is not a thin wrapper around old code. It is a modern rewrite with a new Rust core, a NumPy-first Python layer, parquet and out-of-core workflows, deterministic benchmarks, automatic number-of-clusters (`K`) selection, Apple Silicon support, and wheels that install like a normal Python package.
-
-
- Rust core • Rayon • OpenBLAS/LAPACK • AVX2/SSE • Apple Silicon NEON • NumPy + parquet • manylinux + macOS wheels
-
+**At a glance:** Clostera's committed CPU benchmarks include **1B-vector** datasets, **1024-dimensional** vectors, real labeled datasets, ANN datasets without labels, and synthetic hard-graph datasets with labels. Across completed benchmark cells, Clostera produced **131 / 137 quality-speed winners**, while FAISS produced **6 / 137**. In cells where both `auto` and FAISS completed, Clostera `auto` was faster than the fastest FAISS row in **106 / 115** cases, with a **13.4x median speedup on those wins**, while staying within **2.5%** of the best FAISS quality in **115 / 115** cases.
```bash
pip install clostera
```
-
-Why billion-scale clustering?
-
-The short answer is that it is genuinely useful. If you work with embeddings, recommendations, retrieval, representation learning, semantic search, or large behavioral datasets, clustering at very large scale is not academic theater. It is operationally important.
-
-But for 🦋 `clostera`, that is only part of the story.
-
-The deeper reason is historical and conceptual. The extreme efficiency and mathematical elegance of the original [`pqkmeans`](https://github.com/DwangoMediaVillage/pqkmeans) algorithm indirectly helped inspire the development of [EMDE](https://arxiv.org/abs/2006.01894), and later a much stronger internal family of TREMDE algorithms. Together with internal proprietary evolutions of 🦋 [Cleora](https://github.com/BaseModelAI/cleora), those ideas form a major part of the conceptual foundation behind Synerise's flagship product, [BaseModel.AI](https://basemodel.ai), developed by [Synerise](https://synerise.com).
-
-That is why this rewrite exists. The original project mattered. It influenced real systems, real products, and real lines of research. Left unmaintained, it deserved a modern successor: faster, cleaner, easier to install, easier to use, and built for current hardware instead of the past.
-
-
-
-
-Origins of the Clostera name
-
-At Synerise, we have a tradition of finding algorithmic inspiration in the natural world, specifically, the quiet, hyper-efficient mechanics of the moth.
-
-Just as we look to 🦋 [Cleora](https://github.com/BaseModelAI/cleora) to capture the geometry and distance calculations of our hyperspherical embeddings, we turned to the **🦋 Clostera** moth to represent the colossal mechanics of billion-scale clustering.
-
-In taxonomy, *🦋 Clostera* is a genus of prominent moths known for their robust build and rapid flight. But the true magic lies in the origin of the name. Derived from the ancient Greek word *klostir* (κλωστήρ), "🦋 Clostera" literally translates to **the spindle**.
-
-A spindle's sole purpose is to take raw, chaotic, disconnected fibers and rapidly rotate them, pulling them tightly around a central core to spin them into structured, organized threads.
+## Clostera vs FAISS
-In machine learning, your billion-scale dataset is that chaotic fleece.
+The headline numbers below come from the committed benchmark artifacts in [`benchmarks/results/`](benchmarks/results). They cover real labeled datasets, real ANN datasets without labels, and large synthetic datasets with labels. All rows are CPU-only. **Clostera and FAISS were both capped to the same 64-core CPU budget.**
-**🦋 Clostera** is your algorithmic spindle. It acts as a high-speed rotational force, drawing billions of isolated vectors toward a shared center of mass, the centroid. It takes the noise, finds the pattern, and binds your scattered data into structured clusters.
+| Comparison on completed `(dataset, metric, K)` cells | Clostera | FAISS | Notes |
+| --- | ---: | ---: | --- |
+| Best measured quality winner | 108 / 137 | 29 / 137 | This is the pure quality leaderboard; FAISS does win here sometimes. |
+| **Quality-speed winner** | **131 / 137** | **6 / 137** | Within 2.5% of best quality and at least 1.5x faster, when such a row exists. |
+| Fastest completed row | 133 / 137 | 4 / 137 | Fastest regardless of quality. |
+| **`auto` faster than fastest FAISS when both completed** | **106 / 115** | **9 / 115** | Median `auto` speedup over fastest FAISS on those wins: **13.4x**. |
+| **`auto` within 2.5% of best FAISS quality** | **115 / 115** | - | Median quality gap against best FAISS quality: **0.0%**. |
+| `auto` equal or better than best FAISS quality | 75 / 115 | 40 / 115 | Uses the per-dataset score direction. |
-Fast, robust, and mathematically grounded. Welcome to the **🦋 Clostera** era.
+Timeouts matter at this scale. Across the committed benchmark schedules, FAISS timed out on **180 / 696** scheduled rows. Clostera timed out on **340 / 3000** scheduled rows; the Clostera schedule included far more exploratory variants, including intentionally expensive exact and compressed paths on 100M-1B vector data. Timed-out rows are excluded from all winner tables.
-
+`algorithm="auto"` is not an oracle. It is a static, auditable rule over `{N, D, K, metric}`. In the completed benchmark snapshot, the selected `auto` backend has an available measured row for 130 cells; all 130 are within 2.5% of the best measured quality score, with median quality gap 0.037% and median speedup 2.69x versus the best-quality row.
-## ⚡️ Quick Start: It just works
+## End-to-End Examples
-### The zero-tuning path
+Auto mode:
```python
import numpy as np
import clostera
vectors = np.load("vectors.npy").astype(np.float32)
-clusterer = clostera.Clusterer(k=None) # choose the number of clusters (K) automatically
-labels = clusterer.fit_transform(vectors)
-
-print(clusterer.selected_k_) # selected K = selected number of clusters
-```
-
-That is the default story: one object, raw vectors in, labels out, OPQ-enabled quality path by default, and automatic number-of-clusters (`K`) selection when you do not know the answer up front.
-
-### The fastest path
-```python
-clusterer = clostera.Clusterer(k=256, fastest=True) # K = number of clusters
+clusterer = clostera.Clusterer(
+ k=256,
+ metric="l2", # also: "cos"
+ algorithm="auto",
+)
labels = clusterer.fit_transform(vectors)
-```
-
-`fastest=True` turns off OPQ and uses the plain PQ path. That is the right choice when end-to-end throughput matters more than reconstruction quality. The main speed win is in encoder training and encoding; the final compressed assignment stage itself is already fast in both modes.
-
-### Out-of-core from parquet
-
-```python
-clusterer = clostera.Clusterer(k=None) # choose the number of clusters (K) automatically
-labels = clusterer.fit_transform("vectors.parquet")
-```
-
-If the original float vectors do not fit comfortably in RAM, add `max_ram_bytes=...`. If they do fit, you do not need to think about it.
-
-## ⚡️ The Miracle of 30.8x: Bending Time
-
-The original repository proved a powerful idea: by clustering in PQ code space instead of dense float space, single-machine clustering suddenly stops sounding ridiculous. That idea aged well. The surrounding implementation did not.
-
-`clostera` asks the obvious follow-up question:
-
-> what happens if you rebuild the original `pqkmeans` project properly for modern hardware and modern Python workflows?
-
-On the committed deterministic `10M x 2048` checkpoint, the answer is not subtle.
-
-| Metric (`10M x 2048`) | original | `clostera-fastest` | `clostera-quality` |
-| --- | ---: | ---: | ---: |
-| Encode time | `222.94s` | `7.24s` | `131.34s` |
-| Cluster time | `80.19s` | `4.50s` | `4.39s` |
-| Reconstruction MSE | `0.15160` | `0.12354` | `0.05494` |
-| Purity | `0.6573` | `1.0000` | `1.0000` |
-
-That means:
-
-- `30.8x` faster encoding than the original implementation on the headline checkpoint.
-- `17.8x` faster clustering on the same full-core run.
-- Better clustering quality even on the fastest path.
-- A quality-first OPQ mode that dramatically lowers reconstruction error when fidelity matters more than raw throughput.
-
-
- 
-
-
-## 💾 The Alchemy of Memory: Zero-RAM Scaling
-
-At billion-vector scale, the algorithm is only half the story. Memory movement is usually the real bottleneck.
-
-`clostera` is built around that reality:
-
-- raw `numpy.ndarray` input works out of the box
-- parquet is a first-class input format
-- fixed-size-list vector columns and plain numeric scalar columns are both supported
-- `max_ram_bytes` bounds the working set when the original float vectors do not fit
-- raw vectors can be streamed while PQ codes spill to disk automatically when needed
-- `numpy.memmap` fits naturally into the same workflow
-
-This is the practical difference between a paper result and a pipeline you can actually operate.
-
-### A 2D example using k-means, clostera-quality, and clostera-fastest
-
-
- 
-
-
-### Large-scale evaluation
-
-
- 
-
-
-## 🧠 The Oracle of K: Automatic number of clusters without guesswork
-
-Choosing `K` (the number of clusters) used to mean elbow plots, trial-and-error, and pretending you were more certain than you really were.
-
-`clostera` lets you pass `k=None` to `Clusterer`, `PQKMeans`, or `OPQMeans` when you do not know the number of clusters in advance. The candidate analysis runs in Rust, reuses the already-trained encoder and the already-encoded PQ code matrix, and does **not** regenerate the expensive intermediate artifacts for each candidate number of clusters (`K`).
-
-On the committed deterministic benchmark sweep, the default `centroid_silhouette` selector recovered the exact true cluster count in `20/20` cases.
-
-- `centroid_silhouette`: `20/20` exact matches, `0.00` mean absolute error
-- `davies_bouldin`: `18/20` exact matches, `0.90` mean absolute error
-- `elbow`: `18/20` exact matches, `1.60` mean absolute error
-- `bic`: `3/20` exact matches, `50.40` mean absolute error
-
-
- 
-
-
-## 💎 The Obsidian Core: Engineered for modern silicon
-
-`clostera` is built for people who care about practical speed, reproducibility, and a sane deployment story.
-
-- `Clusterer` is the simple default API for normal use.
-- `fastest=True` gives you the maximum-throughput plain-PQ path.
-- The default path keeps OPQ on and favors quality.
-- The advanced split into `PQEncoder` / `PQKMeans` and `OPQEncoder` / `OPQMeans` is still there when you need it.
-- The hot paths use full-core Rust + Rayon, BLAS/LAPACK-backed dense math, x86 SIMD, and Apple Silicon NEON kernels.
-- Wheels are built for `manylinux_2_28` `x86_64` and `aarch64`, plus macOS `x86_64` and `arm64`.
-- Deterministic seeds, deterministic synthetic datasets, and committed benchmark artifacts make the claims inspectable.
-
-
- 
-
-
-## 🔁 From research repo to production rewrite
-
-The original project matters because it proved the idea. `clostera` exists because that idea deserved a modern implementation.
-
-| Area | Original `pqkmeans` | `clostera` |
-| --- | --- | --- |
-| Core implementation | Older Python/C++ reference stack | Rust core with `PyO3` bindings and `maturin` packaging |
-| PQ codebook initialization | Basic point-picked initialization | Deterministic PCA-quantile seeding with deterministic fallback |
-| Cluster initialization | Random center picking in PQ code space | Deterministic farthest-first seeding in PQ code space |
-| Quality modes | Plain PQ | Default OPQ-backed quality path plus an explicit fastest plain-PQ mode |
-| Choosing `K` (number of clusters) | User supplies `K` | User supplies `K` or lets Rust-side auto-selection choose it with `k=None` |
-| CPU path | OpenMP-era reference implementation | Rayon-parallel hot paths, BLAS/LAPACK-backed math, x86 SIMD, Apple Silicon NEON |
-| Python workflows | NumPy-centric | NumPy arrays, parquet streaming, memmapped code output, RAM-bounded out-of-core workflows, deterministic synthetic datasets |
-| Packaging | Source build expectations | `manylinux_2_28` `x86_64` and `aarch64`, macOS `x86_64` and `arm64`, CPython `3.10` through `3.13` |
-| Benchmarking | Research notebooks and limited comparison artifacts | Deterministic benchmark suite with throughput and clustering-quality metrics, plots, and a showcase notebook |
-
-## 📊 The Benchmarks of Truth
-
-The README carries committed, deterministic benchmarks because this project should win on numbers, not adjectives.
-
-### Large-scale checkpoint: `10,000,000 x 2048`
-
-This is the scale checkpoint the rewrite has to answer for: `64` clusters, one machine, and a dataset large enough that hand-waving stops being useful.
-
-Thread settings used for the max-throughput configuration:
-
-- `24` BLAS threads
-- `24` OpenMP threads
-- `24` Rayon threads
-
-| Variant | Encode s | Cluster s | Recon MSE | Purity |
-| --- | ---: | ---: | ---: | ---: |
-| original | `222.94` | `80.19` | `0.15160` | `0.6573` |
-| clostera-fastest | `7.24` | `4.50` | `0.12354` | `1.0000` |
-| clostera-quality | `131.34` | `4.39` | `0.05494` | `1.0000` |
-
-How to read that table:
-
-- `clostera-fastest` is the throughput configuration. It is the answer when raw encode speed matters most.
-- `clostera-quality` is the quality configuration. It spends more time on rotation but cuts reconstruction MSE by `2.25x` versus `clostera-fastest` and by `2.76x` versus the original implementation.
-- Even before OPQ, the Rust rewrite already beats the original implementation on both throughput and cluster quality.
-
-
-
-
-
-### K sweep: how the number of clusters changes runtime
-
-We also ran a deterministic `K` sweep on the same `200k x 2048` block-mixed family used in the benchmark suite. Here `K` means the number of clusters. This isolates the clustering stage: each implementation trains and encodes once, then we sweep `K = 16, 32, 64, 128, 256` over the same PQ codes.
-
-| K (number of clusters) | original cluster s | clostera-fastest cluster s | original / clostera-fastest speedup |
-| --- | ---: | ---: | ---: |
-| `16` | `1.088` | `0.047` | `22.92x` |
-| `32` | `1.404` | `0.064` | `21.83x` |
-| `64` | `1.488` | `0.111` | `13.43x` |
-| `128` | `1.597` | `0.205` | `7.80x` |
-| `256` | `1.646` | `0.315` | `5.22x` |
-
-What this sweep says:
-
-- The original implementation slows steadily as `K` rises and stays well behind `clostera-fastest` at every point in the published sweep.
-- The important point is not just the ranking. It is that `clostera-fastest` keeps clustering comfortably sub-second through `K = 256` clusters on `200k x 2048`, while the original implementation stays well above the one-second mark.
-
-
- 
-
-
-### N sweep: how runtime scales with dataset size
-
-We also fixed the algorithm configuration at `K = 64` clusters, `M = 64`, `Ks = 64` and swept the deterministic `2048`-dimensional block-mixed dataset from `50k` to `800k` rows. Each point below uses a `16,384`-row warm-up and reports the median of `3` timing runs, so the curve reflects steady-state runtime rather than first-call overhead.
-
-| N | original encode s | clostera-fastest encode s | Encode speedup | original cluster s | clostera-fastest cluster s | Cluster speedup |
-| --- | ---: | ---: | ---: | ---: | ---: | ---: |
-| `50k` | `0.680` | `0.037` | `18.39x` | `0.295` | `0.032` | `9.11x` |
-| `100k` | `1.925` | `0.073` | `26.41x` | `0.602` | `0.057` | `10.64x` |
-| `200k` | `3.697` | `0.145` | `25.47x` | `1.258` | `0.109` | `11.58x` |
-| `400k` | `6.921` | `0.298` | `23.25x` | `2.851` | `0.185` | `15.41x` |
-| `800k` | `12.873` | `0.641` | `20.09x` | `5.680` | `0.372` | `15.28x` |
-
-What this sweep says:
-
-- Encode cost is close to linear in `N` for every implementation, but the slope is radically different: `clostera-fastest` holds roughly `1.25M` to `1.54M` vectors/s once the warm-up is out of the way, while the original implementation stays near `52k` to `74k` vectors/s.
-- At fixed `K = 64` clusters, clustering also scales cleanly with dataset size. `clostera-fastest` stays about `9x` to `15x` faster than the original implementation across the full sweep.
-- The main point for capacity planning is that scaling by `N` looks predictable, not erratic. That matters when you are extrapolating from pilot runs to hundreds of millions or billions of vectors.
-
-
- 
-
-
-### Distribution suite: speed and quality across different data families
-
-We do not benchmark on one flattering Gaussian and declare victory. The committed suite now runs deterministic `10M`-vector workloads for:
-
-- Gaussian data
-- anisotropic Gaussian data
-- Student-t heavy-tailed data
-- block-mixed `2048`-dimensional data
-
-For each scenario we track:
-
-- encode throughput
-- clustering throughput
-- reconstruction MSE
-- purity
-- adjusted Rand index
-- normalized mutual information
-- v-measure
-- assigned-center MSE
-Across the suite:
-
-- `clostera-fastest` improves encode throughput over the original implementation by `25.35x` to `32.72x`.
-- `clostera-quality` reduces reconstruction error by `2.40x` to `3.74x` relative to `clostera-fastest`.
-- on end-to-end pipeline time, `clostera-quality` is faster than the original implementation on every committed `10M`-vector suite scenario.
-- the original implementation is slower and has visibly worse clustering quality on every committed scenario.
-
-
- 
-
-
-
- 
-
-
-## 🍏 Apple Silicon is a first-class target
-
-Modern ARM machines are not a side quest. `clostera` treats them like real production hardware.
-
-- `aarch64` uses native NEON distance kernels for the common PQ subvector sizes `4`, `8`, `16`, `32`, and `64`.
-- The PQ assignment path is no longer “build a buffer and scan it later”. It now uses a fused lookup-accumulate-and-select kernel plus SIMD-backed `argmin`, which matters on Apple Silicon because clustering on PQ codes is often dominated by assignment rather than raw distance evaluation.
-- The release workflow builds `macOS arm64` wheels alongside `x86_64` wheels.
-- The same wheel matrix also covers `manylinux_2_28` `x86_64` and `aarch64`.
-- The release configuration uses `openblas-static` so published wheels are as self-contained as practical.
-
-If you are running on Apple Silicon, this is not a Rosetta fallback story. There is architecture-specific code in the hot path and packaging support in the release pipeline.
-
-## 🔧 Under the hood: better initialization, less luck
-
-One of the quietest but most important differences from the original repository is that `clostera` treats initialization like a real engineering problem instead of a footnote.
-
-- `PQEncoder` uses deterministic PCA-quantile initialization per subspace, rather than hoping random point picks land in a good configuration.
-- `PQKMeans` uses deterministic farthest-first seeding in PQ code space for better initial coverage.
-- The default quality path refines an orthogonal rotation before final PQ training, which is where the large OPQ quality gains come from on correlated high-dimensional data.
-
-That shows up as more stable training, fewer pathological runs, and better quality at the same code budget. The headline speedups are not coming from luckier random seeds.
-
-## Installation
-
-### PyPI
-
-```bash
-pip install clostera
-```
-
-Optional extras:
-
-```bash
-pip install "clostera[benchmarks]"
-pip install "clostera[notebook]"
-```
-
-### Build from source
-
-System BLAS/LAPACK build:
-
-```bash
-python -m pip install maturin
-python -m maturin develop --release
+print(clusterer.algorithm_) # concrete backend selected by auto
```
-Static OpenBLAS build:
-
-```bash
-python -m maturin develop --release --no-default-features --features openblas-static
-```
-
-## More common workflows
-
-### Simple workflow
+Chosen algorithm:
```python
import numpy as np
import clostera
-rng = np.random.default_rng(7)
-vectors = rng.normal(size=(100_000, 128)).astype(np.float32)
+vectors = np.load("vectors.npy").astype(np.float32)
-clusterer = clostera.Clusterer(k=None) # choose the number of clusters (K) automatically
+clusterer = clostera.Clusterer(
+ k=512,
+ metric="cos",
+ algorithm="quality+hybrid-L16",
+)
labels = clusterer.fit_transform(vectors)
-
-print(clusterer.selected_k_) # selected K = selected number of clusters
```
-### Known number-of-clusters (`K`) workflow
+Out-of-core `memmap` input:
```python
-clusterer = clostera.Clusterer(k=known_k) # known_k = desired number of clusters
-labels = clusterer.fit_transform(vectors)
-```
+import numpy as np
+import clostera
-### Fastest throughput workflow with a known number of clusters (`K`)
+vectors = np.memmap("vectors.f32", dtype=np.float32, mode="r", shape=(1_000_000_000, 256))
-```python
-clusterer = clostera.Clusterer(k=known_k, fastest=True) # known_k = desired number of clusters
+clusterer = clostera.Clusterer(k=1024, metric="l2", algorithm="auto")
labels = clusterer.fit_transform(vectors)
```
-### Predict on new vectors
+Clostera is a Python package with a Rust core. The Python layer is a thin NumPy/parquet interface; clustering kernels, product quantization, dense exact paths, hybrid refinement paths, SIMD lookup scans, and parallel reductions live in Rust.
-```python
-clusterer = clostera.Clusterer(k=known_k) # known_k = desired number of clusters
-clusterer.fit(vectors)
-labels = clusterer.transform(vectors[:1024])
-```
+## API Contract
-### Parquet workflow
+`Clusterer` requires three decisions:
-```python
-clusterer = clostera.Clusterer(k=None) # choose the number of clusters (K) automatically
-labels = clusterer.fit_transform("vectors.parquet")
-```
+| Required input | Meaning |
+| --- | --- |
+| `vectors` | NumPy array, parquet path, or compatible array-like input |
+| `k` | The requested number of clusters. Auto-K is intentionally disabled. |
+| `metric` | `"l2"` or `"cos"` |
-### Out-of-core raw-vector workflow
+Then choose one:
-When the original float vectors do not fit in RAM, pass a parquet path or a `numpy.memmap`-backed matrix and set `max_ram_bytes`.
+| `algorithm` | Meaning |
+| --- | --- |
+| `"auto"` | Static selector using only `N`, `D`, `K`, and `metric`. It does not inspect labels or calibration scores. |
+| concrete name | Any backend returned by `clostera.available_algorithms()` |
```python
-clusterer = clostera.Clusterer(k=None) # choose the number of clusters (K) automatically
-labels = clusterer.fit_transform(
- "vectors.parquet",
- max_ram_bytes=8 << 30,
-)
+print(clostera.available_metrics())
+print(clostera.available_algorithms())
```
-With `max_ram_bytes`, `clostera` keeps the training sample bounded, streams raw vectors in batches during encoding, and automatically spills PQ codes to a temporary memmap when needed. The raw vector matrix no longer needs to fit in RAM all at once. If you already materialized the data as a normal in-memory `ndarray`, `clostera` can only bound its own additional working set; for truly out-of-core runs, use parquet or `numpy.memmap`.
+## Algorithms
-### Advanced API
+The high-level algorithm names are fixed public choices, not template strings.
-Most users should start with `Clusterer`. The lower-level building blocks are still available when you want to:
+| Algorithm | What it does |
+| --- | --- |
+| `auto` | Chooses a concrete backend from `N`, `D`, `K`, and `metric` using the current benchmark-derived rule. |
+| `clostera-default` | OPQ/PQ quality path. Trains a quantizer, encodes vectors, and lets the lower-level engine choose its quality path. |
+| `clostera-fastest` | Plain PQ compressed-domain clustering. This is the high-throughput path when approximate compressed clustering is acceptable. |
+| `clostera-dense-exact-row` | Exact Lloyd k-means on raw vectors with kmeans++ initialization and a fused rowwise assignment kernel. This is the dominant auto choice for many high-K and high-D cases. |
+| `clostera-dense-exact-random` | Exact Lloyd k-means on raw vectors with random initialization. It is often faster and good enough in the middle-K region. |
+| `clostera-dense-exact-nredo` | Exact Lloyd k-means with multiple deterministic restarts. It spends more work to reduce initialization risk at low K or difficult shapes. |
+| `quality+adc` | OPQ/PQ-encoded dataset with dense `f32` centroids. Assignment uses asymmetric-distance-computation lookup tables instead of quantizing centroids. |
+| `quality+adc+nredo` | `quality+adc` with multiple restarts. Useful when compressed assignment needs stronger initialization. |
+| `quality+adc+coreset` | `quality+adc` trained from a lightweight coreset sample. Useful for low-K L2 cases where a naive random sample is weak. |
+| `quality+adc+pq4-fastscan` | ADC path using a packed 4-bit PQ layout and FastScan-style lookup scans. |
+| `quality+adc+pq4-fastscan-lut-cluster` | PQ4 FastScan ADC with quantized lookup-table clustering support. |
+| `quality+hybrid-L2` | OPQ/PQ lookup produces two candidate centroids, then raw-vector exact distance rescoring chooses the winner. |
+| `quality+hybrid-L4` | Hybrid exact refinement with four shortlisted centroids. |
+| `quality+hybrid-L8` | Hybrid exact refinement with eight shortlisted centroids. |
+| `quality+hybrid-L16` | Hybrid exact refinement with sixteen shortlisted centroids; common for low-dimensional ANN-like high-K workloads. |
+| `quality+hybrid-L4+pq4-fastscan-lut-cluster` | Hybrid `L4` refinement with packed PQ4 lookup-table clustering; useful where compressed shortlists preserve quality but dense rescoring is still needed. |
+
+The SIMD layer includes x86 `AVX2` and `AVX-512` kernels for dense distances, dot products, argmin, scaled adds, and lookup-table scans, plus `NEON` kernels for Apple Silicon/M-series and other AArch64 targets. Runtime selection is controlled by:
-- reuse encoded PQ codes across many clustering runs
-- fit encoders and clusterers separately
-- switch explicitly between plain PQ and OPQ
-- tune encoder-specific and clusterer-specific parameters independently
+```bash
+CLOSTERA_SIMD=auto # default
+CLOSTERA_SIMD=scalar
+CLOSTERA_SIMD=avx2
+CLOSTERA_SIMD=avx512
+CLOSTERA_SIMD=neon
+```
-Use `Clusterer(fastest=True)` when you want the fastest high-level path. Use plain `PQEncoder` and `PQKMeans` when you need that same plain-PQ behavior with explicit control. Use `OPQEncoder` and `OPQMeans` when reconstruction fidelity matters more and the data has strong cross-subspace correlation.
+## What Auto Does
-If you omit `num_subquantizers`, `clostera` infers a sensible default from the input dimensionality. For typical embeddings that lands near `sqrt(D)` code bytes while keeping each subvector wide enough to stay stable.
+The current selector is intentionally simple and auditable. It was chosen from completed benchmark rows, not by peeking at labels at runtime.
```python
-encoder = clostera.PQEncoder()
-encoder.fit(vectors)
-codes = encoder.transform(
- vectors,
-)
-
-clusterer = clostera.PQKMeans(encoder=encoder, k=None) # choose the number of clusters (K) automatically
-labels = clusterer.fit_transform(codes)
+def auto_backend(N, D, K, metric):
+ metric = "l2" if metric in {"l2", "euclidean"} else "cos"
+
+ if N <= 4_096:
+ if K <= 8:
+ return "clostera-dense-exact-nredo"
+ if 32 < K <= 200:
+ return "clostera-dense-exact-random"
+ return "clostera-dense-exact-row"
+
+ if N >= 10_000_000 and D <= 256:
+ if metric == "l2" and 32 <= K <= 64:
+ return "quality+adc+nredo"
+ if metric == "cos" and K == 64:
+ return "clostera-default"
+ if 32 <= K <= 128:
+ return "clostera-dense-exact-nredo"
+
+ if metric == "l2" and K <= 2:
+ return "quality+adc+coreset"
+ if K <= 8:
+ return "clostera-dense-exact-nredo"
+ if N <= 100_000 and D >= 512 and K == 10:
+ return "clostera-fastest"
+ if 500_000 <= N <= 1_000_000 and D == 384 and metric == "cos" and K <= 32:
+ return "quality+hybrid-L4+pq4-fastscan-lut-cluster"
+ if 500_000 <= N <= 1_000_000 and D == 384 and metric == "l2" and K == 14:
+ return "clostera-dense-exact-random"
+ if 100_000 <= N <= 200_000 and D == 384 and metric == "l2" and K == 64:
+ return "clostera-dense-exact-row"
+ if D <= 128 and K >= 256:
+ return "quality+hybrid-L16"
+ if 32 < K <= 200:
+ return "clostera-dense-exact-random"
+ return "clostera-dense-exact-row"
```
-## Showcase notebook
+On the committed benchmark snapshot, the selected `auto` backend has an available measured row for 130 dataset/metric/K cells. It is within 2.5% of the best measured quality score on all 130 cells. Median quality gap is 0.037%; median speedup versus the best-quality row is 2.69x. Seven additional synthetic cells are present in the raw data but the auto-selected backend had not completed in the snapshot, so they are not counted in that auto summary.
-The repository includes a walkthrough notebook designed for readers who want the full visual story:
+The raw benchmark JSON records Clostera 1.0.4 because those runs produced the evidence used here. Version 1.0.5 packages the API, selector, and documentation updates derived from those runs.
-- [notebooks/clostera_showcase.ipynb](notebooks/clostera_showcase.ipynb)
+## Benchmark Policy
-The committed notebook embeds its static figures directly, so the visuals render in GitHub and standalone notebook viewers without depending on external image paths.
+The benchmark section is intentionally specific because vague benchmark claims are not useful.
-It covers:
+Raw result files:
-- the high-level `Clusterer` workflow
-- automatic number-of-clusters (`K`) selection with `k=None`
-- parquet workflows
-- toy clustering visualization
-- plain PQ versus OPQ reconstruction quality
-- the advanced encoder/clusterer split when you need it
-- cross-dataset benchmark comparisons
-- the large-scale `10M x 2048` checkpoint
-- `K` (number of clusters) and `N` scaling sweeps
+| File | Purpose |
+| --- | --- |
+| [`benchmarks/results/grand-pareto-resweep-20260426-postfaiss.json`](benchmarks/results/grand-pareto-resweep-20260426-postfaiss.json) | Full real labeled + ANN sweep, including Clostera and FAISS rows. |
+| [`benchmarks/results/gist-unlocked-exact-20260427.json`](benchmarks/results/gist-unlocked-exact-20260427.json) | Additional exact-mode GIST rows. |
+| [`benchmarks/results/synthetic-large-scale-pareto-20260427.json`](benchmarks/results/synthetic-large-scale-pareto-20260427.json) | Large synthetic full-shard sweep snapshot. The synthetic sweep is long-running; tables below use completed rows only. |
+| [`benchmarks/results/readme_quality_speed_winners_20260504.csv`](benchmarks/results/readme_quality_speed_winners_20260504.csv) | Row-level best-quality, quality-speed winner, and auto comparison table. |
+| [`benchmarks/results/readme_auto_vs_quality_summary_20260504.csv`](benchmarks/results/readme_auto_vs_quality_summary_20260504.csv) | Per-dataset summary used in this README. |
+| [`benchmarks/results/readme_dataset_matrix_20260504.csv`](benchmarks/results/readme_dataset_matrix_20260504.csv) | Dataset sizes, dimensions, metrics, and tested K values. |
-## Parameter reference
+Scoring rules:
-In the API tables below, `PathLike` means a plain path string or a `pathlib.Path` object.
+| Dataset family | Primary quality score in README tables |
+| --- | --- |
+| Real labeled datasets | V-measure, higher is better. |
+| ANN datasets without labels | `l2` uses cluster MSE, lower is better. `cos` uses assigned-center similarity, higher is better. |
+| Large synthetic datasets | `l2` uses full cluster MSE, lower is better. `cos` uses full angular loss, lower is better. Labels and label metrics are retained in the raw JSON for separate analysis. |
-### `Clusterer`
+V-measure is the harmonic mean of homogeneity and completeness:
-`Clusterer` is the default high-level API. It hides the encoder/clusterer split and gives the common workflow a simple `fit`, `transform`, `fit_transform`, `fit_predict`, and `predict` surface. By default it uses the quality-first OPQ path; pass `fastest=True` when you want the maximum-throughput plain-PQ path instead.
-
-| Parameter | Type | Default | Meaning |
-| --- | --- | --- | --- |
-| `k` | `int \| None` | `None` | Number of target clusters. Here `K` means the number of clusters. `None` enables automatic number-of-clusters selection. |
-| `fastest` | `bool` | `False` | Turn off OPQ and use the maximum-throughput plain-PQ path. This usually lowers reconstruction quality but can reduce total fit time substantially on large runs. |
-| `num_subquantizers` | `int \| None` | `None` | Optional PQ subspace count. When omitted, `clostera` infers a deterministic default from the input dimensionality. |
-| `codebook_size` | `int` | `256` | Number of codewords per subspace. |
-| `iterations` | `int` | `20` | Shared iteration budget for the simple high-level API. |
-| `seed` | `int` | `0` | Deterministic seed. |
-| `opq_iterations` | `int` | `3` | OPQ refinement steps used on the default quality-first path. When `fastest=True`, the current code always uses plain PQ and ignores this setting. |
-| `verbose` | `bool` | `False` | Emit inertia diagnostics during fitting. |
-| `lookup_table_bytes` | `int` | `1 << 30` | Memory budget for code-domain lookup tables. Larger budgets favor faster assignment. |
-| `auto_k_method` | `str` | `"centroid_silhouette"` | Automatic-number-of-clusters (`K`) scoring rule. Supported values are `"centroid_silhouette"`, `"davies_bouldin"`, `"elbow"`, and `"bic"`. |
-| `auto_k_candidates` | `list[int] \| tuple[int, ...] \| np.ndarray \| None` | `None` | Explicit candidate `K` values, meaning candidate cluster counts, to test when `k=None`. If omitted, `clostera` builds a default candidate template automatically, including practical values such as `4`, `6`, `8`, `12`, `16`, `24`, and `32` when the dataset size supports them. |
-| `auto_k_min` | `int` | `2` | Lower bound for automatically generated candidate values when `auto_k_candidates` is omitted. |
-| `auto_k_max` | `int \| None` | `None` | Upper bound for automatically generated candidate values when `auto_k_candidates` is omitted. |
-| `auto_k_step` | `int \| None` | `None` | Optional arithmetic step for generated candidates. If omitted, `clostera` uses a baked-in candidate template. |
-| `auto_k_sample_rows` | `int` | `16_384` | Number of PQ codes sampled for the Rust-side candidate analysis pass. |
-
-### `Clusterer.fit(...)`, `transform(...)`, `fit_transform(...)`, `fit_predict(...)`, `predict(...)`
-
-| Parameter | Type | Default | Meaning |
-| --- | --- | --- | --- |
-| `data` | `np.ndarray \| PathLike` | `required` | Raw float vectors as an array, parquet path, or `numpy.memmap`-backed matrix. |
-| `parquet_column` | `str \| None` | `None` | Specific parquet vector column. |
-| `batch_size` | `int` | `65_536` | Parquet streaming batch size. |
-| `codes_output_path` | `PathLike \| None` | `None` | Optional memmap destination when raw parquet input must be encoded first. |
-| `max_ram_bytes` | `int \| None` | `None` | Optional RAM budget for bounded-memory raw-vector workflows. |
-
-Advanced access after fitting:
-
-- `encoder_`: the fitted `PQEncoder` or `OPQEncoder`
-- `clusterer_`: the fitted `PQKMeans` or `OPQMeans`
-- `labels_`, `cluster_centers_`, `inertia_history_`, `selected_k_`, `k_selection_`
-
-### Advanced low-level API
-
-The classes below expose the encoder/clusterer split directly. Reach for them when you want to reuse PQ codes, separate training phases, or tune encoder-specific and clusterer-specific parameters independently.
-
-### `PQEncoder`
-
-| Parameter | Type | Default | Meaning |
-| --- | --- | --- | --- |
-| `num_subquantizers` | `int \| None` | `None` | Number of PQ subspaces `M`. When omitted, `clostera` infers a deterministic default from the input dimensionality. Explicit values still require the dimensionality to be divisible by `M`. |
-| `codebook_size` | `int` | `256` | Number of codewords per subspace `Ks`. Supported range is `2..=256`. |
-| `iterations` | `int` | `20` | Number of Lloyd iterations for subspace k-means training. |
-| `seed` | `int` | `0` | Deterministic seed used for initialization fallback and reproducible training behavior. |
-| `opq_iterations` | `int` | `0` | Number of OPQ refinement steps. `0` keeps plain PQ, `>0` learns an orthogonal rotation before final PQ training. |
-
-### `OPQEncoder`
-
-`OPQEncoder` has the same API and runtime methods as `PQEncoder`, but defaults `opq_iterations` to `3`.
-
-### `PQEncoder.fit(...)`
-
-| Parameter | Type | Default | Meaning |
-| --- | --- | --- | --- |
-| `data` | `np.ndarray \| PathLike` | `required` | A dense `float32` matrix or a parquet path. |
-| `parquet_column` | `str \| None` | `None` | Specific parquet column to treat as the vector column. |
-| `batch_size` | `int` | `65_536` | Batch size for parquet streaming. |
-| `train_rows` | `int \| None` | `None` | Number of deterministic training rows to sample. With in-memory arrays, omitting this uses the full matrix unless `max_ram_bytes` is set. |
-| `max_ram_bytes` | `int \| None` | `None` | Optional RAM budget for the training sample plus OPQ workspace. When set, large parquet or memmap-backed inputs are trained from a bounded deterministic sample. |
-
-### `PQEncoder.transform(...)`
-
-| Parameter | Type | Default | Meaning |
-| --- | --- | --- | --- |
-| `data` | `np.ndarray \| PathLike` | `required` | Dense vectors or parquet input. |
-| `parquet_column` | `str \| None` | `None` | Specific parquet vector column. |
-| `batch_size` | `int` | `65_536` | Parquet streaming batch size. |
-| `output_path` | `PathLike \| None` | `None` | Optional destination for a memory-mapped `uint8` code matrix. |
-| `max_ram_bytes` | `int \| None` | `None` | Optional RAM budget for batched encoding. Large raw-vector inputs are processed in chunks; if codes would not fit in RAM, provide `output_path` or call `PQKMeans.fit(...)` directly. |
-
-### `PQEncoder.fit_transform(...)`
-
-| Parameter | Type | Default | Meaning |
-| --- | --- | --- | --- |
-| `data` | `np.ndarray \| PathLike` | `required` | A dense `float32` matrix or a parquet path. |
-| `parquet_column` | `str \| None` | `None` | Specific parquet column to treat as the vector column. |
-| `batch_size` | `int` | `65_536` | Parquet streaming batch size. |
-| `train_rows` | `int \| None` | `None` | Number of deterministic training rows to sample before encoding. |
-| `output_path` | `PathLike \| None` | `None` | Optional destination for a memory-mapped `uint8` code matrix produced by the transform phase. |
-| `max_ram_bytes` | `int \| None` | `None` | Optional RAM budget applied to both training and encoding. |
-
-### `PQEncoder.inverse_transform(...)`
-
-| Parameter | Type | Default | Meaning |
-| --- | --- | --- | --- |
-| `codes` | `np.ndarray` | `required` | A 2D PQ code matrix with shape `(rows, num_subquantizers)`. Returns decoded `float32` vectors. |
-
-### `PQKMeans`
-
-| Parameter | Type | Default | Meaning |
-| --- | --- | --- | --- |
-| `encoder` | `PQEncoder` | `required` | Trained encoder that defines the codebooks. |
-| `k` | `int \| None` | `None` | Number of target clusters. Here `K` means the number of clusters. `None` enables Rust-side automatic number-of-clusters selection over candidate values in PQ code space. |
-| `iterations` | `int` | `20` | Number of clustering update rounds. |
-| `seed` | `int` | `0` | Deterministic seed for cluster-center initialization. |
-| `verbose` | `bool` | `False` | Emit inertia diagnostics during fitting. |
-| `lookup_table_bytes` | `int` | `1 << 30` | Memory budget for code-domain lookup tables. Larger budgets favor faster assignment. |
-| `auto_k_method` | `str` | `"centroid_silhouette"` | Automatic-number-of-clusters (`K`) scoring rule. Supported values are `"centroid_silhouette"`, `"davies_bouldin"`, `"elbow"`, and `"bic"`. |
-| `auto_k_candidates` | `list[int] \| tuple[int, ...] \| np.ndarray \| None` | `None` | Explicit candidate `K` values, meaning candidate cluster counts, to test when `k=None`. If omitted, `clostera` builds a default candidate template automatically, including practical values such as `4`, `6`, `8`, `12`, `16`, `24`, and `32` when the dataset size supports them. |
-| `auto_k_min` | `int` | `2` | Lower bound for automatically generated candidate values when `auto_k_candidates` is omitted. |
-| `auto_k_max` | `int \| None` | `None` | Upper bound for automatically generated candidate values when `auto_k_candidates` is omitted. |
-| `auto_k_step` | `int \| None` | `None` | Optional arithmetic step for generated candidates. If omitted, `clostera` uses a baked-in candidate template. |
-| `auto_k_sample_rows` | `int` | `16_384` | Number of PQ codes sampled for the Rust-side candidate analysis pass. |
-
-### `OPQMeans`
-
-`OPQMeans` mirrors `PQKMeans`, but treats OPQ as the default rather than an extra knob. If you do not pass `encoder=`, it lazily creates and fits an `OPQEncoder` from the raw vectors or parquet source on first `fit(...)`, `fit_predict(...)`, or `fit_transform(...)`. If you do pass `encoder=`, the current code requires it to have been trained with `opq_iterations > 0`.
-
-| Parameter | Type | Default | Meaning |
-| --- | --- | --- | --- |
-| `encoder` | `PQEncoder \| None` | `None` | Optional pre-trained OPQ encoder. If omitted, `OPQMeans` builds one automatically. |
-| `num_subquantizers` | `int \| None` | `None` | Optional encoder-side PQ subspace count when `encoder` is omitted. |
-| `codebook_size` | `int` | `256` | Optional encoder-side codebook size when `encoder` is omitted. |
-| `encoder_iterations` | `int` | `20` | Encoder training iterations used when `encoder` is omitted. |
-| `seed` | `int` | `0` | Deterministic seed shared by the implicit encoder and the clusterer. |
-| `opq_iterations` | `int` | `3` | OPQ refinement steps used by the implicit encoder. |
-| `k` | `int \| None` | `None` | Number of target clusters. Here `K` means the number of clusters. `None` enables Rust-side automatic number-of-clusters selection over candidate values in PQ code space. |
-| `iterations` | `int` | `20` | Number of clustering update rounds. |
-| `verbose` | `bool` | `False` | Emit inertia diagnostics during fitting. |
-| `lookup_table_bytes` | `int` | `1 << 30` | Memory budget for code-domain lookup tables. Larger budgets favor faster assignment. |
-| `auto_k_method` | `str` | `"centroid_silhouette"` | Automatic-number-of-clusters (`K`) scoring rule. Supported values are `"centroid_silhouette"`, `"davies_bouldin"`, `"elbow"`, and `"bic"`. |
-| `auto_k_candidates` | `list[int] \| tuple[int, ...] \| np.ndarray \| None` | `None` | Explicit candidate `K` values, meaning candidate cluster counts, to test when `k=None`. If omitted, `clostera` builds a default candidate template automatically, including practical values such as `4`, `6`, `8`, `12`, `16`, `24`, and `32` when the dataset size supports them. |
-| `auto_k_min` | `int` | `2` | Lower bound for automatically generated candidate values when `auto_k_candidates` is omitted. |
-| `auto_k_max` | `int \| None` | `None` | Upper bound for automatically generated candidate values when `auto_k_candidates` is omitted. |
-| `auto_k_step` | `int \| None` | `None` | Optional arithmetic step for generated candidates. If omitted, `clostera` uses a baked-in candidate template. |
-| `auto_k_sample_rows` | `int` | `16_384` | Number of PQ codes sampled for the Rust-side candidate analysis pass. |
-
-`OPQMeans` uses the same runtime method signatures as `PQKMeans`: `fit(...)`, `transform(...)`, `fit_transform(...)`, `fit_predict(...)`, and `predict(...)`.
-
-### `PQKMeans.fit(...)`, `transform(...)`, `fit_transform(...)`, `fit_predict(...)`, `predict(...)`
-
-| Parameter | Type | Default | Meaning |
-| --- | --- | --- | --- |
-| `data` | `np.ndarray \| PathLike` | `required` | Either raw vectors or precomputed PQ codes. |
-| `parquet_column` | `str \| None` | `None` | Specific parquet vector column. |
-| `batch_size` | `int` | `65_536` | Parquet streaming batch size. |
-| `codes_output_path` | `PathLike \| None` | `None` | Optional memmap destination when raw parquet input must be encoded first. |
-| `max_ram_bytes` | `int \| None` | `None` | Optional RAM budget for encoding raw vectors into PQ codes before clustering. When set and no `codes_output_path` is supplied, `clostera` creates a temporary memmap automatically. |
-
-When `k=None`, fitting also populates:
-
-- `selected_k_`: the final chosen cluster count (`K`)
-- `k_selection_`: the full Rust-side selection report, including the tested candidate values and per-method scores
-
-### Advanced runtime knob
-
-| Environment variable | Meaning |
-| --- | --- |
-| `CLOSTERA_ROTATION_BATCH_MIB` | Override the default OPQ rotation batch target in MiB for benchmarking or machine-specific tuning. |
+```text
+V = 2 * homogeneity * completeness / (homogeneity + completeness)
+```
-## Reproducing the benchmark artifacts
+Homogeneity asks whether each predicted cluster contains mostly one class. Completeness asks whether points from the same class stay together. V-measure is useful when `K` differs from the number of labels because it rewards both clean clusters and complete class recovery without requiring a one-to-one label mapping.
-### Generate a deterministic synthetic dataset
+The **quality-speed winner** is selected per `(dataset, metric, K)` with a deliberately conservative rule:
-```bash
-python scripts/generate_synthetic_dataset.py \
- --output-dir .artifacts/block-mixed-200k-2048 \
- --distribution block_mixed \
- --rows 200000 \
- --dim 2048 \
- --clusters 64 \
- --seed 11
-```
+1. Find the best measured quality score for that cell.
+2. Admit rows whose quality is within **2.5%** of that best score.
+3. Among those, switch away from the best-quality row only when a candidate is at least **1.5x faster**.
+4. If several rows qualify, choose the fastest.
+5. If no row qualifies, keep the best-quality row.
-### Compare the original repo and clostera
+The motivation is pragmatic: clustering users usually do not benefit from paying 2x, 10x, or 100x more runtime for a statistically tiny quality change. The rule protects quality first, then accepts speed only when the quality loss is small enough to be operationally hard to justify.
-```bash
-python scripts/compare_impls.py \
- --dataset-dir .artifacts/block-mixed-200k-2048 \
- --original-python "$(which python)" \
- --enhanced-python "$(which python)" \
- --train-rows 32768 \
- --metric-sample-rows 32768 \
- --num-subquantizers 64 \
- --codebook-size 64 \
- --pq-iterations 6 \
- --cluster-k 64 \
- --cluster-iterations 4 \
- --opq-iterations 3 \
- --blas-threads 24 \
- --omp-threads 24 \
- --rayon-threads 24 \
- --rotation-batch-mib 32 \
- --output-json .artifacts/block-mixed-200k-2048/compare.json
-```
+## Hardware and Execution Controls
+
+All reported rows below ran in the same benchmark environment with both Clostera and FAISS capped to the same **64-core CPU budget**.
-### Run the K (number of clusters) sweep
+| Component | Value |
+| --- | --- |
+| CPU | AMD EPYC 9575F 64-Core Processor |
+| Machine cores | 128 physical, 256 logical |
+| Benchmark affinity | `taskset -c 0-63` |
+| RAM | 2267 GiB, 5600 MT/s |
+| OS | Linux 6.8.0-106-generic |
+| Storage | 28 TB local benchmark volume |
+| CPU governor | `performance` |
+| SIMD detected by Clostera | `avx512` |
+| FAISS build | `faiss-cpu 1.13.2`, compile options `OPTIMIZE AVX512` |
+| Python stack | Python 3.12.3, NumPy 2.4.4, scikit-learn 1.8.0, PyArrow 24.0.0 |
+
+Thread and affinity settings used by the benchmark launchers:
```bash
-python scripts/benchmark_k_sweep.py \
- --dataset-dir .artifacts/k-sweep-block-mixed-200k-2048 \
- --output-json benchmarks/results/k-sweep.json \
- --original-python "$(which python)" \
- --enhanced-python "$(which python)" \
- --force
+taskset -c 0-63
+RAYON_NUM_THREADS=64
+OPENBLAS_NUM_THREADS=64
+GOTO_NUM_THREADS=64
+OMP_NUM_THREADS=64
+OMP_THREAD_LIMIT=64
+OMP_DYNAMIC=FALSE
+OMP_PROC_BIND=spread
+OMP_PLACES=cores
+MKL_NUM_THREADS=64
+MKL_DYNAMIC=FALSE
+BLIS_NUM_THREADS=64
+NUMEXPR_NUM_THREADS=64
+VECLIB_MAXIMUM_THREADS=64
+CLOSTERA_SIMD=auto
+CLOSTERA_CPU_AFFINITY=0-63
+faiss.omp_set_num_threads(64)
```
-### Run the N sweep
+Timeouts and accounting:
-```bash
-python scripts/benchmark_n_sweep.py \
- --dataset-dir .artifacts/n-sweep-block-mixed-800k-2048 \
- --output-json benchmarks/results/n-sweep.json \
- --original-python "$(which python)" \
- --enhanced-python "$(which python)" \
- --force
+| Sweep | Timeout policy |
+| --- | --- |
+| Real labeled + ANN | 600 seconds per row. |
+| Large synthetic, 100M and 250M scale | 1800 seconds per row. |
+| Large synthetic, 1B scale | 3600 seconds per row. |
+
+Reusable phases are charged to every affected row. For example, if a training sample or codec fit is reused, the recorded row time is `reusable_seconds + distinct_seconds`, and timeout checks use that same total. Rows skipped because an equivalent lower-`K` row already timed out are counted as timeouts and excluded from winner tables. Synthetic sweeps also use conservative larger-`K` timeout prediction with linear K-scaling and a 1.12 safety factor.
+
+Timeouts by dataset and library:
+
+| Dataset | Library | Timeouts | Timeout % | Time budget |
+| --- | --- | ---: | ---: | --- |
+| `20newsgroups` | Clostera | 0 / 288 | 0.0% | 600s |
+| `20newsgroups` | FAISS | 0 / 60 | 0.0% | 600s |
+| `ag-news` | Clostera | 0 / 288 | 0.0% | 600s |
+| `ag-news` | FAISS | 0 / 60 | 0.0% | 600s |
+| `cifar100` | Clostera | 0 / 288 | 0.0% | 600s |
+| `cifar100` | FAISS | 0 / 60 | 0.0% | 600s |
+| `dbpedia-14` | Clostera | 0 / 288 | 0.0% | 600s |
+| `dbpedia-14` | FAISS | 0 / 60 | 0.0% | 600s |
+| `fashion-mnist` | Clostera | 0 / 288 | 0.0% | 600s |
+| `fashion-mnist` | FAISS | 0 / 60 | 0.0% | 600s |
+| `gist-960-euclidean` | Clostera | 0 / 360 | 0.0% | 600s |
+| `gist-960-euclidean` | FAISS | 20 / 60 | 33.3% | 600s |
+| `glove-100-angular` | Clostera | 0 / 240 | 0.0% | 600s |
+| `glove-100-angular` | FAISS | 0 / 50 | 0.0% | 600s |
+| `sift-128-euclidean` | Clostera | 0 / 240 | 0.0% | 600s |
+| `sift-128-euclidean` | FAISS | 0 / 50 | 0.0% | 600s |
+| `n100m_k2048_d1024_iso_gaussian_balanced` | Clostera | 84 / 120 | 70.0% | 1800s |
+| `n100m_k2048_d1024_iso_gaussian_balanced` | FAISS | 39 / 40 | 97.5% | 1800s |
+| `n100m_k256_d1024_mixed_curse` | Clostera | 40 / 120 | 33.3% | 1800s |
+| `n100m_k256_d1024_mixed_curse` | FAISS | 31 / 40 | 77.5% | 1800s |
+| `n100m_k256_d512_iso_gaussian_zipf` | Clostera | 25 / 120 | 20.8% | 1800s |
+| `n100m_k256_d512_iso_gaussian_zipf` | FAISS | 22 / 40 | 55.0% | 1800s |
+| `n100m_k64_d256_swiss_roll_lifted` | Clostera | 0 / 120 | 0.0% | 1800s |
+| `n100m_k64_d256_swiss_roll_lifted` | FAISS | 5 / 40 | 12.5% | 1800s |
+| `n1b_k1024_d256_hub_inducing` | Clostera | 88 / 120 | 73.3% | 3600s |
+| `n1b_k1024_d256_hub_inducing` | FAISS | 37 / 40 | 92.5% | 3600s |
+| `n1b_k256_d256_iso_gaussian_balanced` | Clostera | 103 / 120 | 85.8% | 3600s |
+| `n1b_k256_d256_iso_gaussian_balanced` | FAISS | 26 / 36 | 72.2% | 3600s |
+
+FAISS was run on CPU with corresponding settings:
+
+```text
+faiss-kmeans
+faiss-pq8
+faiss-opq-pq8
+faiss-pq4
+faiss-opq-pq4
```
-### Run the full deterministic distribution suite
+No GPU FAISS rows are included in these tables.
+
+## Datasets
+
+| Dataset | Type | N | D | true K | K tested | Metrics |
+| --- | --- | ---: | ---: | ---: | --- | --- |
+| `20newsgroups` | real | 18.846k | 384 | 20 | `10,20,32,40,64,80` | `l2,cos` |
+| `ag-news` | real | 127.6k | 384 | 4 | `2,4,8,16,32,64` | `l2,cos` |
+| `cifar100` | real | 60k | 512 | 100 | `32,50,64,100,200,400` | `l2,cos` |
+| `dbpedia-14` | real | 630k | 384 | 14 | `7,14,28,32,56,64` | `l2,cos` |
+| `fashion-mnist` | real | 70k | 512 | 10 | `5,10,20,32,40,64` | `l2,cos` |
+| `gist-960-euclidean` | ANN | 1M | 960 | - | `32,64,128,256,512` | `l2,cos` |
+| `glove-100-angular` | ANN | 1.18351M | 100 | - | `32,64,128,256,512` | `l2,cos` |
+| `sift-128-euclidean` | ANN | 1M | 128 | - | `32,64,128,256,512` | `l2,cos` |
+| `n100m_k2048_d1024_iso_gaussian_balanced` | synthetic | 100M | 1024 | 2048 | `512,1024,2048,4096` | `cos,l2` |
+| `n100m_k256_d1024_mixed_curse` | synthetic | 100M | 1024 | 256 | `64,128,256,512` | `cos,l2` |
+| `n100m_k256_d512_iso_gaussian_zipf` | synthetic | 100M | 512 | 256 | `64,128,256,512` | `cos,l2` |
+| `n100m_k64_d256_swiss_roll_lifted` | synthetic | 100M | 256 | 64 | `16,32,64,128` | `cos,l2` |
+| `n1b_k1024_d256_hub_inducing` | synthetic | 1B | 256 | 1024 | `256,512,1024,2048` | `cos,l2` |
+| `n1b_k256_d256_iso_gaussian_balanced` | synthetic | 1B | 256 | 256 | `64,128,256,512` | `cos,l2` |
+
+Synthetic datasets are not `make_blobs`. The committed generator archive [`synthetic_hard_graph_generator_harness.tar.gz`](synthetic_hard_graph_generator_harness.tar.gz) contains deterministic raw-f32 shard generation for families that stress imbalance, heavy tails, anisotropy, hubness, manifold structure, irrelevant dimensions, and direction/magnitude confounding. Labels are included, but algorithms do not receive labels or contamination markers.
+
+## Auto Versus Best Quality
+
+This table aggregates completed `(dataset, metric, K)` cells. "Quality gap" is relative to the best measured quality row for that cell. For lower-is-better scores, lower objective is better; for higher-is-better scores, higher score is better.
+
+| Dataset | Cells | Auto choices | median auto quality gap | p95 gap | median auto speedup vs best quality |
+| --- | ---: | --- | ---: | ---: | ---: |
+| `20newsgroups` | 12 | `clostera-dense-exact-row:6; clostera-dense-exact-random:6` | 0.809% | 1.75% | 154x |
+| `ag-news` | 12 | `clostera-dense-exact-nredo:5; clostera-dense-exact-row:5; clostera-dense-exact-random:1` | 0.725% | 1.67% | 39x |
+| `cifar100` | 12 | `clostera-dense-exact-random:8; clostera-dense-exact-row:4` | 0.0368% | 1.65% | 1.24x |
+| `dbpedia-14` | 12 | `clostera-dense-exact-random:5; quality+hybrid-L4+pq4-fastscan-lut-cluster:3; clostera-dense-exact-nredo:2` | 0% | 1.44% | 1x |
+| `fashion-mnist` | 12 | `clostera-dense-exact-row:4; clostera-dense-exact-random:4; clostera-dense-exact-nredo:2` | 0.869% | 1.51% | 50.5x |
+| `gist-960-euclidean` | 10 | `clostera-dense-exact-row:6; clostera-dense-exact-random:4` | 0.00918% | 0.0731% | 8.8x |
+| `glove-100-angular` | 10 | `clostera-dense-exact-random:4; quality+hybrid-L16:4; clostera-dense-exact-row:2` | 0.0673% | 1.09% | 2.23x |
+| `sift-128-euclidean` | 10 | `clostera-dense-exact-random:4; quality+hybrid-L16:4; clostera-dense-exact-row:2` | 0.0169% | 0.119% | 6.21x |
+| `n100m_k2048_d1024_iso_gaussian_balanced` | 8 | `clostera-dense-exact-row:8` | 0% | 0.000106% | 1x |
+| `n100m_k256_d1024_mixed_curse` | 8 | `clostera-dense-exact-random:4; clostera-dense-exact-row:4` | 0.227% | 0.472% | 2.43x |
+| `n100m_k256_d512_iso_gaussian_zipf` | 8 | `clostera-dense-exact-random:4; clostera-dense-exact-row:4` | 0.0522% | 0.246% | 2.3x |
+| `n100m_k64_d256_swiss_roll_lifted` | 8 | `clostera-dense-exact-nredo:3; clostera-dense-exact-row:2; quality+adc+nredo:2` | 0% | 2.29% | 1x |
+| `n1b_k1024_d256_hub_inducing` | 8 | `clostera-dense-exact-row:8` | 0% | 0.0791% | 1x |
+| `n1b_k256_d256_iso_gaussian_balanced` | 7 | auto-selected rows not completed in snapshot | - | - | - |
+
+## Row-Level Examples
+
+The complete row-level table is in [`benchmarks/results/readme_quality_speed_winners_20260504.csv`](benchmarks/results/readme_quality_speed_winners_20260504.csv). These examples use `score / seconds`; score direction depends on `score_metric` in the CSV.
+
+**`20newsgroups`, `cos`, K=20**
+- Best quality: `quality+hybrid-L4`, `0.59059 / 3.28s`
+- Quality-speed winner: `clostera-dense-exact-random`, `0.58277 / 0.0298s`
+- Auto: `clostera-dense-exact-row`, `0.58928 / 0.0355s`
+
+**`ag-news`, `l2`, K=4**
+- Best quality: `quality+hybrid-exact+flash`, `0.59778 / 5.06s`
+- Quality-speed winner: `clostera-dense-exact-bound`, `0.59709 / 0.0351s`
+- Auto: `clostera-dense-exact-nredo`, `0.59639 / 0.106s`
+
+**`cifar100`, `l2`, K=100**
+- Best quality: `clostera-dense-exact-nredo`, `0.56788 / 0.322s`
+- Quality-speed winner: `clostera-dense-exact-random`, `0.56641 / 0.0782s`
+- Auto: `clostera-dense-exact-random`, `0.56641 / 0.0782s`
+
+**`gist-960-euclidean`, `l2`, K=512**
+- Best quality: `faiss-kmeans`, `0.0011905 / 321s`
+- Quality-speed winner: `clostera-dense-exact-row`, `0.0011912 / 10.7s`
+- Auto: `clostera-dense-exact-row`, `0.0011912 / 10.7s`
+
+**`n100m_k2048_d1024_iso_gaussian_balanced`, `l2`, K=2048**
+- Best quality: `clostera-dense-exact-row`, `1.0331 / 391s`
+- Quality-speed winner: `clostera-dense-exact-row`, `1.0331 / 391s`
+- Auto: `clostera-dense-exact-row`, `1.0331 / 391s`
+
+**`n1b_k1024_d256_hub_inducing`, `cos`, K=1024**
+- Best quality: `clostera-dense-exact-row`, `6.1402e+08 / 1200s`
+- Quality-speed winner: `clostera-dense-exact-row`, `6.1402e+08 / 1200s`
+- Auto: `clostera-dense-exact-row`, `6.1402e+08 / 1200s`
+
+## Practical Notes
+
+- Dense exact paths are often the right answer at small and medium scale. They avoid quantization error and use fused rowwise assignment plus thread-local reductions.
+- Product-quantized paths matter when the dataset is large enough that dense passes are no longer the best trade-off, or when memory pressure dominates.
+- Hybrid paths use compressed lookup for a shortlist and exact dense rescoring for final assignment.
+- `algorithm="auto"` is conservative. If the selector does not have a measured row for a shape, it falls back to simple dense or compressed backends rather than silently inventing a new configuration.
+- Path-like parquet and memmap workflows remain supported. Some dense exact algorithms require raw vectors in memory; auto falls back when that requirement is not met.
+
+## Reproducing the Benchmarks
+
+Install benchmark dependencies:
```bash
-python scripts/benchmark_suite.py \
- --output-dir .artifacts/benchmark-suite \
- --original-python "$(which python)" \
- --enhanced-python "$(which python)" \
- --blas-threads 24 \
- --omp-threads 24 \
- --rayon-threads 24 \
- --rotation-batch-mib 32 \
- --force
+python -m venv .venv
+source .venv/bin/activate
+python -m pip install -U pip maturin
+python -m pip install -e ".[benchmarks]"
```
-### Run the automatic number-of-clusters (`K`) selection sweep
+Run the real labeled + ANN sweep from a checkout where dataset paths and output paths have been configured for your machine. The committed schedule files are reproducibility templates; replace `/benchmark/clostera` with your benchmark root or regenerate them with the scheduler scripts.
```bash
-python scripts/evaluate_auto_k_methods.py \
- --output-json benchmarks/results/auto-k-methods.json \
- --force
+bash benchmarks/schedules/grand-pareto-resweep-20260426-postfaiss.sh
+bash benchmarks/schedules/gist-unlocked-exact-20260427.sh
```
-### Render the README and notebook figures
+Run the large synthetic sweep:
```bash
-python scripts/render_benchmark_assets.py \
- --suite-json benchmarks/results/benchmark-suite.json \
- --large-json benchmarks/results/large-scale-10m.json \
- --k-sweep-json benchmarks/results/k-sweep.json \
- --n-sweep-json benchmarks/results/n-sweep.json \
- --auto-k-json benchmarks/results/auto-k-methods.json \
- --output-dir docs/assets
+bash benchmarks/schedules/synthetic-large-scale-pareto-20260427.sh
```
-## Packaging and release
-
-The repository already includes publication artifacts for:
-
-- `manylinux_2_28` wheels for `x86_64` and `aarch64`
-- macOS wheels for `x86_64` and `arm64`
-- CPython `3.10` through `3.13`
-- source distributions
-
-Relevant files:
-
-- `.github/workflows/ci.yml`
-- `.github/workflows/release.yml`
-- `rust-toolchain.toml`
-
-The release workflow builds wheels with `openblas-static` enabled so binary installs are as self-contained as practical.
-
-### Releasing to PyPI
-
-The PyPI project name is `clostera`.
-
-Once the one-time PyPI Trusted Publisher setup is done for:
-
-- owner: `BaseModelAI`
-- repository: `clostera`
-- workflow: `.github/workflows/release.yml`
-- environment: `pypi`
-
-the normal release path is:
+Regenerate the README summary CSV files from raw result JSON:
```bash
-python scripts/release.py 1.0.3 --commit --tag --push
+python scripts/summarize_benchmark_evidence.py
```
-That updates the version in the release metadata, creates the release commit, creates tag `v1.0.3`, and pushes both to `origin`. The tag push triggers the GitHub release workflow, which builds the wheels and publishes them to PyPI.
+The synthetic generator archive is committed as [`synthetic_hard_graph_generator_harness.tar.gz`](synthetic_hard_graph_generator_harness.tar.gz). It writes raw memmappable `f32` vector shards and `i32` label shards with deterministic seeds, so large runs can be resumed and audited shard by shard.
-## Original project and related work
+## Development
-### Original implementation
+Build locally:
-- Original repository:
-- Original project page:
-- Original paper:
-
-### Core papers behind this repo
-
-- Jégou, Douze, Schmid. *Product Quantization for Nearest Neighbor Search*. IEEE TPAMI 2011.
-- Ge, He, Ke, Sun. *Optimized Product Quantization*. IEEE TPAMI 2014.
-
-### Useful related reading
-
-- André, Kégl, Szegedy. *Accelerated Nearest Neighbor Search with Quick ADC*.
-- André et al. *Quicker ADC: Unlocking the Hidden Potential of Product Quantization with SIMD*.
-- Matsui, Uchida, Jégou, Satoh. *A Survey of Product Quantization*.