[Aminsed]

Summary

v0.1.1 returned a single repeated value for every input on any system with
NumPy 2.x. The CUDA kernels were computing correct values; the binding threw
them away. This PR fixes that, replaces the two hot kernels, and rebuilds the
packaging, tests, CI and documentation around them.

Verified on an RTX A6000 (sm_86, CUDA 12.5): 45 Python tests and 11 C++ tests
pass, ruff is clean, and a clean pip install from the tree works.

The headline bug

python/CMakeLists.txt pinned pybind11 v2.10.4, which predates NumPy 2.0's
change to the array descriptor layout. Built against NumPy ≥ 2, every returned
array had a stride of 0, so all elements aliased element 0:

cuLOF   : 1 unique value;   [9.4945 9.4945 9.4945 9.4945 ...]
sklearn : 500 unique values; [9.4945 9.9042 8.9527 10.4107 ...]

The kernels were fine — the value being broadcast is the correct score for
point 0. CMake now requires pybind11 ≥ 2.12, and
test_output_is_not_a_constant guards the symptom.

Other correctness fixes

• Every build targeted sm_52. project(LANGUAGES CUDA) seeds
  CMAKE_CUDA_ARCHITECTURES, so the detection function's if(DEFINED ...)
  guard always returned early — an sm_86 A6000 ran Maxwell code.
• normalize=True was undefined behaviour: cudaMemcpyHostToDevice on a
  device pointer, and cudaHostRegister on device memory.
• Non-contiguous input silently gave wrong answers: the binding checked only
  strides[1] before passing the raw pointer, so X[::2] was read as
  contiguous.
• k > 32 raised cudaErrorInvalidValue; zero-variance features produced NaN
  under normalize=True; --use_fast_math contradicted the documented aim of
  matching scikit-learn; BUILD_TESTS=ON could not configure.

Performance

Distances are now a single cuBLAS SGEMM via ‖a-b‖² = ‖a‖² + ‖b‖² - 2a·b on
centred data, replacing an elementwise kernel that accumulated in float64 — a
1:64 rate on consumer Ampere.

k-NN is an exact radix select on the raw float bit patterns. LOF needs the
neighbour set and the k-distance, never their order, so selection suffices —
and the k-th key the select produces is the k-distance. No k limit, coalesced
reads, O(n) per pass. The Gram-to-distance conversion is fused into the first
pass and 11-bit digits cover 32 bits in three passes, so the tile is swept four
times rather than six.

The n×n matrix is no longer materialised; tiling keeps memory linear in n.
n = 200,000 now runs where it previously needed 149 GiB.

samples	scikit-learn	cuLOF	speedup
1,000	0.013 s	0.0002 s	63×
20,000	0.717 s	0.0221 s	33×
50,000	2.970 s	0.0954 s	31×
200,000	16.93 s	1.411 s	12×

RTX A6000 vs 24-thread i9-12900K, k=20, best of 7, context init excluded.

API

LOF is now a genuine drop-in for sklearn.neighbors.LocalOutlierFactor —
n_neighbors, contamination, fit_predict returning ±1 labels,
score_samples, negative_outlier_factor_, offset_. A test asserts label-for-label
equality with scikit-learn. culof.lof(X, k) returns the LOF value itself.

The previous set_k/set_threshold surface is gone, as is the pure-Python
LOF class that __init__.py defined and then immediately overwrote with the
C++ one — so the documented API never existed at runtime.

Accuracy

Agreement is measured, not asserted: median relative difference 4.5×10⁻⁷,
p99 2.4×10⁻⁶, identical top-1% rankings everywhere. The tail is points whose
k-th and (k+1)-th neighbours are closer than float32 can resolve, where the
choice of k-th neighbour is genuinely ambiguous; 0.27% of points, and it never
changed which points were flagged. test_near_ties_are_the_only_disagreement
pins this down rather than hiding it behind a loose tolerance.

Everything else

• Packaging on pyproject.toml + scikit-build-core; setup.py gone; version has one source.
• Python package moved under src/ so tests exercise the installed artefact, not the source tree.
• Removed the duplicated python/ tree, which shipped a second copy of the tests and a top-level module named python.
• Deterministic output: the gather uses a block-wide prefix sum, not atomics, so runs are bit-identical.
• CI across CUDA 11.8/12.4 and Python 3.9–3.12.
• Figures are generated from the benchmark's own JSON. The old README and two figures carried three contradictory speedup numbers (6.75×, 4.6×, 1.57×).

Note

The .github/workflows/ci.yml file needs a token with the workflow scope to
push over HTTPS; this branch was pushed over SSH.