Metadata-Version: 2.4
Name: anydist
Version: 1.0.2
Summary: HNSW approximate nearest-neighbor search for arbitrary data and dissimilarity functions
Author-email: Matteo Dell'Amico <della@linux.it>
License-Expression: BSD-3-Clause
Project-URL: Homepage, https://gitlab.com/bobtables/anydist
Keywords: hnsw,nearest-neighbor,ann,similarity-search,non-metric,custom-distance
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Developers
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Cython
Classifier: Topic :: Scientific/Engineering :: Artificial Intelligence
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Dynamic: license-file

# anydist

Approximate nearest-neighbor search on **arbitrary data with arbitrary
dissimilarity functions**, using Hierarchical Navigable Small World graphs
(HNSW; [Malkov & Yashunin](https://arxiv.org/abs/1603.09320)).

Unlike most HNSW libraries, which index numeric vectors under a fixed set of
metrics, this one accepts any Python objects and any dissimilarity function
you write — including non-metric ones. It also has the usual fast path: a set
of built-in metrics over numeric vectors, computed in SIMD C kernels with no
Python on the hot path, with parallel bulk insertion and parallel batch
queries. Between the two there is a third door: hand it a **C function
pointer** — from numba, ctypes, Cython or an extension module of your own —
and *your* dissimilarity runs on that same GIL-free hot path (see [C
metrics](#c-metrics-numba-ctypes-an-extension-module)).

**The index stays live.** You can add to it and delete from it after it is
built, interleaved with queries, without a rebuild — which not every library
offers: nmslib and annoy are batch indexes, where a point added after the
build costs a full rebuild to become findable (see
[Living indexes](#living-indexes)).

Two implementations:

- **`HNSWFast`** — the one you want: a Cython implementation over typed,
  array-backed data structures. `HNSW` is an alias for it, and a normal install
  builds it.
- **`HNSWPure`** — the pure-Python, dict-based *reference* implementation, no
  dependencies beyond the standard library. It exists to document the
  algorithm readably and to check the fast path against, not as a deployment
  option: with a built-in metric it is ~30x slower to build and ~40x slower to
  query (measured on 5k 20-dimensional vectors), and even with a Python
  callback — where the callback itself dominates — it is ~1.5x slower.

**`HNSW` never falls back to `HNSWPure`.** Without the compiled extension
there is no `HNSW`, and asking for it raises `ImportError` saying so — an
order of magnitude is not something to lose silently to a build that did not
happen. If you want the pure implementation, name it:

```python
from anydist import HNSWPure as HNSW      # deliberately, penalty and all
```

In a hurry? [**Quickstart**](#quickstart) is two code blocks, and the
[**API reference**](#api) documents every constructor argument and method in
one place — go straight there if you would rather read signatures than prose.

This code has been developed together with [fishdbc](https://gitlab.com/bobtables/fishdbc), an incremental density-based clustering algorithm that leverages this library for scalability and allowing arbitrary dissimilarity functions.

## Installation

```sh
pip install anydist
```

or, from a checkout of this repository:

```sh
pip install .
```

Windows is not supported: MSVC, the compiler Windows Pythons are built with,
lacks the atomic builtins the parallel code needs. Under WSL, the Linux wheels
work as usual.

PyPI has prebuilt wheels for the common platforms. Anywhere else, and from a
checkout, pip compiles the Cython extension, so you need a C compiler with
OpenMP; on macOS, `brew install libomp` first (or point `LIBOMP_PREFIX` at
another libomp). That extension is what gives you `HNSWFast`. A failed build
shows up the first time you ask for `HNSW`, since there is no fallback to hide
it, but you can check directly:

```python
import anydist
assert anydist.HNSWFast is not None
```

On x86-64 Linux, the portable build already selects AVX2 kernels at load time
on CPUs that have them, so it is not leaving the obvious speedup on the table.
For a further boost tied to the machine you build on (at the cost of a
non-portable binary):

```sh
HNSW_NATIVE=1 pip install .
```

On macOS the package carries its own copy of the OpenMP runtime, `libomp`, and
a process must never load two: libomp aborts when it finds another copy of
itself already running. An extension of your own that uses OpenMP alongside
anydist should therefore build against this copy rather than bring one:
`anydist.get_openmp_dir()` returns its directory, and its docstring has the
compiler and linker flags. fishdbc does exactly this.

## Quickstart

Arbitrary objects, your own dissimilarity function:

```python
from anydist import HNSW

def jaccard(a, b):
    return 1 - len(a & b) / len(a | b)

index = HNSW(jaccard)                       # the dissimilarity is the 1st argument
index.update(my_sets)                       # insert them all
neighbors = index.search(query_set, k=10)   # [(id, distance), ...], nearest first
```

Ids are insertion order: the `id` in a result is the position of that element
in the sequence you inserted. The index is a read-only sequence of those
elements, so `index[id]` gives one back, `len(index)` counts them, and
iteration and slicing work.

Numeric vectors and a built-in metric — no Python callback on the hot path:

```python
index = HNSW('euclidean')            # vector length inferred from the first row
index.update(X)                      # X: (n, dim) array, inserted in parallel
I, D = index.search_batch(Q, k=10)   # (n_queries, k) ids and distances, in parallel
```

With a built-in metric, `update` inserts across all cores; with a callback
metric it stays serial by default, since only you can say whether your
function is safe to run on several threads — pass `num_threads=N` to opt in.
`add(elem)` is still there for genuinely incremental use, where you interleave
inserts and queries.

`index[id]` works in both modes, but they store elements differently: a
callback metric keeps your objects, while a built-in metric keeps rows of a
typed matrix and hands you back a copy of one — rounded to float32, since that
is the default storage there (see `float32`).

## Multithreading

`update` and `search_batch` both take `num_threads`: `None` (the default)
picks for you, a positive value is a thread count, and **`-1` (or any value
≤ 0) means one thread per CPU**.

The rule is that **anything we wrote defaults to every core, anything you wrote
defaults to one**, because only you can say whether your dissimilarity function
`d` — the callable you passed as `metric` — is thread-safe. What you get once
you opt in depends on the metric, because the GIL does:

- **Built-in metric** — ours, and stateless, so the kernels are `nogil` and the
  threads are real on any interpreter build. Measured on 12 CPUs: 5.5x build /
  5.0x query at 50k x 32, and 4.3x / 3.7x at Fashion-MNIST's 60k x 784. Don't
  read those as an efficiency: the machine's cores are not all the same speed,
  so what a given box can offer is its own question.
- **C metric** (a function pointer you supplied) — the call is `nogil` too, so
  opting in is the only way a *user-written* distance gets real threads on a
  stock interpreter: measured 4.3x at one-per-core (n=4k, dim=16, a
  deliberately slow `jit_metric`). It still defaults to serial, because
  `nogil` buys the *ability* to run your `d` in parallel, not the safety of
  doing so — mutable state reached through `user_data`, a memo table, a
  counter, all race, with no GIL here to accidentally serialise them.
- **Callback metric** — the callback runs under a thread state re-attached per
  distance. On a free-threaded interpreter that is real parallelism; on a
  stock GIL build it pays off only when each distance holds the GIL briefly
  and releases it for a long time — measured 3.0x at 2 threads for a
  40k-dimensional NumPy distance, but *slower than serial* at 128 dimensions,
  where contention and the thread-state attach cost more than the op. Measure
  before believing.

These figures are indicative, not a spec — see [About these
numbers](#about-these-numbers) for the method and its resolution limits.

```python
index.update(X)                              # built-in metric: all cores
index.update(objs)                           # your metric: one thread
index.update(objs, num_threads=-1)           # opt in: "my d is thread-safe"
index.update(objs, num_threads=4)            # opt in, four threads
I, D = index.search_batch(Q, num_threads=1)  # force a single-threaded query
```

`search_batch` is worth calling for a callback metric even serially: it is the
same graph walk as a `search` loop, but it hands back the whole batch as two
`(n_queries, k)` arrays — neighbor ids and their distances — instead of a list
of lists.

`chunksize`, the constructor's or a call's own, interacts with this: workers take chunks of
*consecutive* ids, which is what keeps parallel insertion correct when the
input is order-local (sorted, time-series, grouped records) — two near
neighbors inserted at the same instant never see each other, and the edge
between them is lost. The default adapts to the batch; raise it if your input
is strongly ordered.

### Writing a callback metric that scales

On a free-threaded interpreter there is one rule, and it is not about how you
write the function — it is about what the function *reads*:

> **A parallel metric must not read a shared mutable object on every call.**
> Reading `self.something` or a closed-over variable costs 3-4x at eight
> threads. Read parameters from module globals or default arguments instead.

CPython gives some objects immortal or deferred reference counts, so touching
them from many threads is free: `None`, small ints, literals, module-level
functions, and **all classes** (including ones built at runtime by `type()`,
`namedtuple`, or a `class` statement inside a function). Everything else is
ordinary, and an object that every thread increfs on every distance call turns
one cache line into a contention point. Measured on the same laptop (n=10 000
tuples of 8 floats, `ef=25`, 8 threads, free-threaded 3.13):

| metric written as | speedup at 8 threads |
|---|---|
| module-level `def`, reads only its arguments | **3.53x** |
| closure that **reads** its captured variable | 0.82x — *slower than serial* |
| closure that ignores it | 3.08x |
| bound method that **reads** `self` | 1.19x |
| bound method that ignores `self` | 2.94x |

Note the pairs: same callable form, same arithmetic, only the read differs. The
form is irrelevant — `functools.partial` and lambdas are fine. So a metric with
parameters should carry them like this:

```python
# slow in parallel: `self.p` and the closure cell are one shared object each
class Metric:
    def __call__(self, a, b):
        return sum((x - y) ** self.p for x, y in zip(a, b))

# all of these scale (3.3-3.5x measured):
P = 2
def metric(a, b):                    # module-level global
    return sum((x - y) ** P for x, y in zip(a, b))

def make_metric(p):
    def metric(a, b, p=p):           # default-argument binding
        return sum((x - y) ** p for x, y in zip(a, b))
    return metric
```

Per-element objects are *not* a problem: your data points are thousands of
distinct objects, so the traffic spreads. It is the single shared one that hurts.

**If threading makes things *slower*, suspect numpy.** There is a bug in numpy
([numpy#32298](https://github.com/numpy/numpy/issues/32298)) where an array
that has been through `pickle`, `copy.deepcopy`, or a `joblib` cache ends up
sharing one small internal object that every thread has to synchronise on. The
array is otherwise identical, and single-threaded work is unaffected — but on a
free-threaded build the contention can turn a speedup into a slowdown.

Usually you need not care: `update()` detects and repairs it whenever you hand
it an array. The case it cannot reach is a list of rows you extracted yourself,
because then the library never sees the array. Repair it at the source:

```python
from anydist import canonicalize_dtype
rows = list(canonicalize_dtype(X))       # X came from pickle / joblib / a cache
index.update(rows, num_threads=-1)
```

So: parallel inserts slower than serial, on data that came out of a cache, is
worth one try of the two lines above — see [`BUGS.md`](BUGS.md), which also
lists the other issue worth knowing about.

## Built-in metrics

Passing a name instead of a callable computes the distance in a SIMD C kernel
with no Python on the hot path. Available names (aliases in parentheses):

| Name | Distance | Notes |
|---|---|---|
| `euclidean` (`l2`) | `sqrt(sum (x-y)^2)` | |
| `l2sqr` (`sqeuclidean`) | `sum (x-y)^2` | ranks identically to `l2`, skips the sqrt |
| `cosine` (`cosinesimil`) | `1 - cos_sim` | normalizes for you; the value scipy and sklearn report |
| `angular` (`angulardist`) | `arccos(cos_sim) / pi` | normalizes for you; ranks as `cosine` does, and obeys the triangle inequality |
| `dot` | `1 - dot(x, y)` | normalizes **nothing** — cosine distance only if you pre-normalized |
| `negdotprod` | `-dot(x, y)` | the same space as `dot`, shifted by 1 |
| `l1` (`manhattan`) | `sum \|x-y\|` | |
| `linf` (`chebyshev`) | `max \|x-y\|` | |

The names follow nmslib's dense-vector spaces. `dot` and `negdotprod` differ by
an additive constant, and every decision the search makes is a comparison, so
they build identical graphs and return identical orders — choose between them
by the number you want back. `dot` is non-negative on unit-length rows, which
a downstream consumer of the weights may need; `negdotprod` is what code
ported from nmslib expects.

### Cosine distance

Three names get you cosine-like behavior, and they differ in whether they
normalize for you and in what number comes back. All three return the **same
neighbor lists** on unit-length rows, so this is a choice about values and
cost, not about quality.

**Just want cosine distance? Use `cosine`.** It divides by the length of each
of the two vectors itself, so your rows do not have to be unit-length, and it
returns `1 - cos_sim` — the same value scipy, scikit-learn, nmslib
(`cosinesimil`) and hnswlib report.

```python
index = HNSW('cosine')
index.update(X)                          # any lengths; normalizing is on us
index.search(q, k=10)                    # distances are 1 - cos(q, x), in [0, 2]
```

**Need the triangle inequality? Use `angular`.** Same normalization, but it
returns the angle between the vectors scaled to `[0, 1]` instead of
`1 - cos_sim`. Since a smaller angle is exactly a larger cosine, it ranks
candidates identically to `cosine`; unlike `cosine`, it is a true metric on
directions. The extra `arccos` is not free — it is the only libm call in any
built-in kernel, and it costs roughly a third of query time at 50 dimensions —
so reach for `angular` when something downstream needs metric behavior, not by
default.

**Already normalized, and want the cheapest option? Use `dot`.** It normalizes
nothing, so on unit rows `1 - dot(x, y)` *is* cosine distance, computed with a
single accumulator instead of three.

```python
import numpy as np
from anydist import HNSW

X = np.asarray(X, dtype=np.float64)
X /= np.linalg.norm(X, axis=1, keepdims=True)   # unit rows (drop zero rows first)

index = HNSW('dot')
index.update(X)

q = q / np.linalg.norm(q)                # the query too: dot normalizes nothing
index.search(q, k=10)
```

A zero row has no direction, so `cosine` and `angular` both report distance 0
for it rather than a NaN.

## C metrics: numba, ctypes, an extension module

A built-in metric is fast because it never enters Python; a callback is
flexible because it is Python. This section is for when you want both — and it
asks more of you than the two options above, so skip it unless a Python
callback has turned out to be your bottleneck.

Give `HNSW` a **C function pointer** instead of a Python callable, and your
distance runs with no GIL and no Python on the hot path, exactly where a
built-in metric runs. Parallel insertion and parallel queries then become
possible on a stock interpreter — but you have to ask for them:

```python
index.update(X, num_threads=-1)              # -1 = one thread per core
I, D = index.search_batch(Q, k=10, num_threads=-1)
```

**A C metric still defaults to one thread**, like a Python callback and unlike
a built-in metric. `nogil` buys the *ability* to run your distance in
parallel, not the safety of doing so: state reached through `user_data`, a
memo table or a counter all race, and there is no GIL here to serialise them
by accident. Only you can say whether your `d` is thread-safe, so the opt-in
is yours to make.

The signature is fixed, and it is the one scipy's `LowLevelCallable` users will
recognise:

```c
double f(const double *a, const double *b, Py_ssize_t n, void *user_data);
```

### Writing one in Python: `jit_metric`

Nobody wants to spell that signature in numba by hand, so the decorator
generates it from the two-argument form:

```python
from anydist import HNSW, jit_metric

@jit_metric                          # needs numba installed
def manhattan(a, b):                 # a, b: float64 arrays of length dim
    s = 0.0
    for i in range(a.shape[0]):
        s += abs(a[i] - b[i])
    return s

index = HNSW(manhattan)
index.update(X, num_threads=-1)      # opt in: "my d is thread-safe"
I, D = index.search_batch(Q, k=10, num_threads=-1)
```

It returns a `numba.cfunc`, which `HNSW` recognises by its `.address`. The body
must compile in `nopython` mode — loops over the two arrays are the idiom,
arbitrary Python is not; if you need arbitrary Python, pass the plain function
and pay the GIL per distance. `@jit_metric(cache=True)` persists the compiled
kernel across processes, which matters for a benchmark harness that runs one
cell per process; it requires the function to live in an importable file.

What the pointer buys, measured on the [same 12 CPUs](#about-these-numbers)
(stock-GIL CPython 3.14, n=5000, dim=32, m=16, ef=100, k=10) with the *same*
distance written both ways — the callback being the obvious numpy one-liner:

| metric | threads | build | query |
|---|---|---|---|
| Python callback | 1 | 5.34 s | 1.63 ms/q |
| C metric (`jit_metric`) | 1 | 0.20 s | 0.032 ms/q |
| C metric | 12 | 0.044 s | 0.008 ms/q |

The 12-thread row is opt-in (`num_threads=-1`) — a C metric is serial by
default, see [Multithreading](#multithreading) — and it is safe here because
this metric is a pure function of its two arguments.

26x on build and 50x on query before threading, 4.6x and 4.3x more from it. The
serial ratio is the part that moves with your metric: a callback whose body
does real work narrows it, one written in pure Python widens it. The threading
factor is the part a callback cannot have on a stock build at all.

### Pointing it at a library you already have

The pointer is duck-typed, so nothing here is imported unless you use it. Any
of these shapes is accepted:

- **scipy `LowLevelCallable`** — `.function` is a `PyCapsule`, `.user_data`
  carries parameters. scipy already normalises numba, ctypes and Cython into
  it, so this one shape covers all three.
- **numba `cfunc`** — has `.address`.
- **ctypes function pointer** — how an existing `.so` reaches Python.
- anything else exposing `.address` or a `.function` capsule (plus an optional
  `.user_data`), including a `PyCapsule` from your own extension module.

```python
import ctypes

dbl = ctypes.POINTER(ctypes.c_double)
lib = ctypes.CDLL("./libdtw.so")
lib.dtw_metric.restype = ctypes.c_double
lib.dtw_metric.argtypes = [dbl, dbl, ctypes.c_ssize_t, ctypes.c_void_p]

index = HNSW(lib.dtw_metric)         # DTW, in C, across all cores
```

**Parameters go through `user_data`**, which is read from every one of those
shapes: keep them in a struct (or a single `ctypes` cell) that outlives the
index, and expose its address as `.user_data` next to the pointer.
`jit_metric` does not expose `user_data` — a closure over module-level
constants covers the parameterised case with less ceremony, and a hand-built
`cfunc` or a `LowLevelCallable` is there when you want the pointer itself.

**One thing to know before you take this door: getting the signature wrong
*crashes* the process instead of raising an exception.** How much protection
you get depends on what you hand over:

All three forms declare their signature somewhere, and all three are checked
against `double (const double *, const double *, Py_ssize_t, void *)`, raising
a `CSignatureWarning` if they disagree:

- **A capsule** carries it in its name — a `scipy.LowLevelCallable`, or
  anything else whose `.function` is a named PyCapsule. Spelling is not the
  test: `const`, spacing, parameter names and `Py_ssize_t`/`ssize_t`/
  `npy_intp` all pass; the types have to line up.
- **A ctypes function pointer** carries it as `_restype_`/`_argtypes_`, and a
  **numba cfunc** carries the same on its `.ctypes` view. Either pointer
  spelling is accepted for the vectors — `POINTER(c_double)` or a bare
  `c_void_p`. `jit_metric` builds the signature for you, so that path cannot
  disagree in the first place.

What is genuinely unchecked is a **bare address**: an unnamed capsule, or a
pointer whose declaration you never wrote down. There the crash contract stands
in full.

The check **warns rather than refuses**, because it recognises a list of
spellings rather than parsing C, and an unfamiliar typedef must not be able to
block a pointer that is in fact correct. If that is your situation:

```python
warnings.simplefilter("ignore", anydist.CSignatureWarning)
```

and please open an issue with the spelling, so the list can grow.

One limit is worth knowing before you write the metric rather than after:
**your vectors must all be the same length.** `n` is a single dimension,
inferred from the first row and passed to every call. Variable-length data
stays on the callback path, which imposes nothing on what your objects are.

Otherwise a C metric follows the built-in-metric path rather than the callback
one, which the [API reference](#constructor) covers argument by argument. One
of those defaults is worth a second look, because it is the only place that
path assumes something about *your* function rather than ours: `select` is
chosen from the vector length, on the reasoning that a short vector is cheap
to re-measure. That holds for the built-in kernels, whose work is linear in
the length. If yours is expensive at any length — an alignment, a search, a
simulation — pass **`select='edges'`** and it will read stored weights instead
of calling you O(m²) times per row.

## API

### Constructor

```python
HNSW(metric, m=5, ef=200, m0=None, level_mult=None,
     level_strategy='random', keep_pruned=False, select=None, cache=None,
     float32=None, keep_workspace=False, chunksize=None)
```

- **`metric`** — the dissimilarity, in one of three forms (it need not satisfy
  the metric axioms in any of them — the name follows the scikit-learn/hdbscan
  convention, not a metricity requirement): a callable `d(x, y)` over arbitrary
  data; a built-in name (see above) over numeric vectors, whose length is
  inferred from the first one inserted; or a **C function pointer** over such
  vectors (numba `cfunc`, `scipy.LowLevelCallable`, ctypes pointer, capsule),
  which runs GIL-free like a built-in one — see [C
  metrics](#c-metrics-numba-ctypes-an-extension-module). A pointer is
  recognised *before* callability, so an object that is both takes the fast
  door.
- **`m`, `ef`, `m0`, `level_mult`** — the paper's parameters; see
  [Malkov & Yashunin](https://arxiv.org/abs/1603.09320).
- **`level_strategy`** — how insertion levels are assigned. `'random'` is the
  paper's strategy and the best choice for approximate nearest-neighbor search;
  `'balanced'` is a deterministic variant (an element is promoted when its
  neighborhood is full and no neighbor sits above it) that yields a more
  balanced structure.
- **`keep_pruned`** — the paper's optional `keepPrunedConnections` (Alg. 4,
  lines 15–17). Off by default.
- **`select`** — the neighbor-selection rule.
  `'edges'` answers the paper's pruning heuristic (Alg. 4) from edges the graph
  already has, spending zero extra distance evaluations; `'distances'` is the
  heuristic as written, recomputing the candidate-to-result distances;
  `'naive'` just takes the `m` nearest (Alg. 3). `HNSWFast` defaults to
  `'edges'` for a callback metric, and for a built-in one decides at the first
  vector: `'distances'` while `dim <= 4*m` (recomputing a short vector beats
  fetching a stored weight), `'edges'` above that — a C metric follows that
  same first-vector rule. `HNSWPure` always defaults to `'edges'`, since every
  distance there is a Python call.

  That first-vector rule takes vector length as a stand-in for what a distance
  costs, which is true of the built-in kernels. **A C metric that is expensive
  out of proportion to its vector length should be given `select='edges'`
  explicitly** — otherwise short rows get `'distances'`, which evaluates your
  function O(m²) times per row instead of reading stored weights.
- **`cache`** — reuse, within a single insert or query, distances computed at a
  previous level of the same descent. Defaults on for a callback metric (13%
  faster builds, measured), off for a built-in one — and off for a C metric,
  which is cheap for the same reason — where the check costs about what it
  saves.
- **`float32`** — store built-in-metric vectors as float32. Defaults on for a
  built-in metric, as hnswlib and faiss do: it halves both the vector store and
  the bandwidth of a query's dominant read.
  It does *round the data*, so pass `False` when the exact float64
  distances are the point.
  `True` is **rejected** for the other two metric forms rather than quietly
  overridden: a callback metric keeps Python objects and has no typed store to
  narrow, and a C metric reads its rows as `const double *`, so a float32 store
  would be reinterpreted rather than converted. Both leave the vectors as
  float64, and `index[i]` hands back a copy of the stored row either way.
- **`keep_workspace`** — hold the threads' working memory between calls
  instead of allocating it fresh each time. Only affects `update()` and
  `search_batch()` with more than one thread.

  Each worker needs a set of buffers sized by the number of nodes in the
  index, so **every threaded call pays an allocation proportional to the index
  — not to the work you asked for**. One big `update()` never notices. A loop
  of small ones pays it every time, and it grows as the index does. At 100k
  nodes and four threads:

  | | fresh each call | kept |
  |---|---|---|
  | `update()` in chunks of 200 | 4.01 ms per chunk | **2.46 ms** |
  | `search_batch()` of 64 rows | 2.11 ms | **0.37 ms** |

  Off by default because it is a trade, not a free win: the buffers stay
  allocated — `threads × nodes × 32 bytes`, or 44 with `cache=True`, so 370 to
  500 MB at twelve threads and a million nodes — until you call
  `release_workspace()`. Peak
  memory is the same either way; what changes is the floor between calls. Turn
  it on if you build or query in a stream of small batches, leave it off
  otherwise, and call `release_workspace()` when that phase is over.

  Results are identical either way, and the setting survives pickling while
  the buffers themselves do not.
- **`chunksize`** — how many *consecutive* ids each worker takes in a threaded
  `update()` that does not pass its own; see `update` below for what the
  chunk guards. `None`, the default, adapts it to each batch. `HNSWPure`
  accepts and keeps it, with no parallel insertion to schedule.

### Methods

- **`update(X, num_threads=None, ef=None, chunksize=None)`** — insert many
  elements. `X` is a numeric array for a built-in **or C** metric — anything
  array-like will do, and it is converted to C-contiguous float64, so integer
  rows, float32 rows and lists of lists are all accepted and all stored as
  float64 — or any sequence of objects for a callback metric. Parallel by
  default for a built-in metric, serial by default for the two you supply
  (see [Multithreading](#multithreading)). `chunksize` sets how many
  *consecutive* ids each worker takes: the chunk is what keeps parallel
  insertion correct on order-local input (sorted, time-series, grouped), where
  two near neighbors inserted concurrently would never see each other. `None`
  takes the constructor's `chunksize`, and if that is `None` too the chunk
  adapts to the batch; set it, on the index or per call, if you know your
  input is strongly order-local.
- **`add(elem, ef=None)`** — insert one element, returning nothing; its id is
  the previous `len(index)`.
- **`search(q, k=None, ef=None)`** — the (approximate) `k` closest elements to
  `q` as `(id, distance)` pairs, sorted by distance. `k=None` returns all `ef`
  candidates found. Soft-deleted nodes are excluded (see below).
- **`search_batch(Q, k=None, ef=None, num_threads=None)`** — answer many
  queries at once, optionally in parallel. `Q` is `(n_queries, dim)` for a
  built-in or C metric, or a sequence of objects for a callback one; returns
  `(indices, distances)` of shape `(n_queries, k)` — `k` defaulting to `ef` —
  padded with `-1`/`inf` where fewer than `k` neighbors were found. Both
  classes have it, with the same rows, widths and padding, so
  `I, D = index.search_batch(Q, k)` reads the same either way. `HNSWFast`
  returns numpy arrays and can use threads; `HNSWPure` returns lists of lists,
  having no numpy, and loops over `search`, accepting `num_threads` for
  signature parity and ignoring it as its `update` does.
- **`index[i]`, `len(index)`, iteration, slicing** — the index is a read-only
  sequence of the elements you inserted, in id order (negative ids count from
  the end; a built-in or C metric returns copies of the stored rows). Note this
  is *not* what `index[i]` meant in this package's predecessor, where it
  yielded node `i`'s neighbors; `graphs` is where those live now.
- **`mark_deleted(i)` / `unmark_deleted(i)`** — soft deletion (see below);
  `index.deleted_count` reports how many are currently marked.
- **`add_recording(elem, ef=None, harvest_cap=0)`** — insert one element and
  return every distance computed during the insertion as three parallel
  sequences `(i, j, dist)`, where `i` is always the new id — numpy arrays from
  `HNSWFast`, and from `HNSWPure` the standard library's `array.array` (typecodes
  `'i'`, `'i'`, `'d'`), as compact at 16 bytes a triple. A plain `add` records
  nothing and pays no overhead.
- **`update_recording(X, ..., harvest_cap=0)`** — the batch form, returning the
  same three sequences for the whole insertion. It is `add_recording` in a
  loop, and given the same insertion levels the two harvest the same triples.
  Note that `HNSWFast`'s comes back in per-worker order rather than insertion
  order, so sort by `i` if you need it.

  `harvest_cap` bounds what comes back, identically on both: `0` is off, and a
  positive `N` keeps only each inserted node's `N` shortest triples. It does
  **not** change the index — the same distances are computed and the same graph
  is built either way, so it trades recorded detail for the downstream cost of
  consuming it. Selection is by distance alone, which can keep the wrong bridge
  when a nearer cluster fills the quota first; leave it off for streaming
  workloads with concept drift, where insertion order cannot shuffle that away.
  Any cap `>= 1` keeps every node attached, since a node is only ever matched
  against already-inserted ones.
- **`freeze()`** — give up inserting in exchange for the memory only insertion
  needed (the stored edge weights: `n*m0*8` bytes at the base layer, plus
  the smaller per-node blocks above it). There is no way
  back, as with a `frozenset`. Under
  `select='distances'`, which stores no weights to begin with, it only
  installs the guard. Returns `None`: it freezes the index you call it on
  rather than handing back a frozen copy.
- **`release_workspace()`** — free the threads' working memory that
  `keep_workspace=True` holds on to. Safe to call at any time, including when
  nothing is held; the next threaded call simply allocates again. Present on
  both classes, so code holding either can call it without asking which.
- **`graphs`** — the adjacency as a list of `{node: {neighbor: dist}}`, one per
  level, for inspection.

### Living indexes

An index you can keep changing is not universal among ANN libraries. Adding
2,000 vectors one at a time to a 20k-node index, with a k=10 query every 100
inserts (dim=32, m=16, serial):

| library | insert after build | μs/insert | μs/delete |
|---|---|--:|---|
| **this library** | yes | 61 | `mark_deleted`, 0.1 |
| hnswlib | yes | 53 | yes, 0.1 |
| usearch | yes | 73 | yes, 0.6 |
| faiss (`IndexHNSWFlat`) | yes | 129 | **none** |

nmslib and annoy are not in the table because they are batch indexes and say
so — nmslib: "only static data sets are supported"; annoy: "after calling
build, no more items can be added". A point added later is kept, but the
structure that answers queries was built once and does not know about it, so
getting it in means rebuilding the whole index. That is a different operation
from an incremental insert, not a slower one, and there is nothing here to
compare it against.

### Deletion

`mark_deleted(i)` drops node `i` from results; `unmark_deleted(i)` puts it
back, and `deleted_count` says how many are out. Deletion is soft, as in
hnswlib: the node stops appearing in results but stays a routing waypoint, so
the graph never gets disconnected and recall for everything else is untouched.
The cost is that deleted nodes keep their memory and still cost distance
computations while routing, so an index that is mostly deleted wants a rebuild.

```python
index.mark_deleted(7)
ids, dists = index.search_batch(Q, k=10)   # 7 will not appear
```

The exclusion happens *inside* the search, not over its output: deleted nodes
are still traversed, so what comes back is the true nearest **live**
neighbors — not the unfiltered answer with rows struck out (measured >0.99
recall against brute force over the live set, with half the index deleted).
You get `min(k, ef, live)` results: the search fills its pool with `ef` *live*
candidates before it is willing to stop, so a short list means there really are
that few live nodes, not that the search gave up. (`k > ef` truncates to `ef`,
deletions or not — that is the pool size, so keep `ef` comfortably above `k` as
usual.)

### Persistence

`anydist.save` and `anydist.load` take a path or a binary file object:

```python
import anydist

anydist.save(index, "index.pkl")
index = anydist.load("index.pkl")          # same graph, same results
```

There is no `saves()`/`loads()` pair returning bytes — pass an `io.BytesIO()`,
as with `numpy.save`. Plain `pickle` (or `dill`) still works on an index too;
`save`/`load` exist for the one thing pickle cannot do on its own.

**The metric is the only part that may not travel.** A built-in name or a
module-level function pickles; a lambda, a closure or a C function pointer does
not, and an address would be meaningless in the loading process anyway.
Everything expensive — the vectors, the adjacency, the deletions — always
serialises, so rebuilding would be a poor answer. Instead, save without the
metric and give it back on load:

```python
anydist.save(index, "index.pkl", include_metric=False)
index = anydist.load("index.pkl", metric=my_cfunc)
```

Including the metric is the default **and raises if it cannot be done**, so
`save` never quietly writes a file that `load(f)` cannot open. You get a
self-contained file or an exception naming the way out — never a surprise in
the process that loads it.

**Speed.** Saving and loading are close to memory-bandwidth bound, so an
index round-trips far faster than it builds. Measured on the machine described
under [About these numbers](#about-these-numbers), writing to a real file in
a temporary directory — as every library in the comparison does, so the
figures are like-for-like:

| index | on disk | save | load |
|---|--:|--:|--:|
| Fashion-MNIST, 60k x 784 | 213 MB | 0.18 s | 0.17 s |
| Gaussian, 100k x 128 | 92 MB | 0.09 s | 0.06 s |
| Gaussian, 50k x 32 | 14 MB | 0.02 s | 0.01 s |

Roughly 1.2 GB/s each way — the same order as a memory copy, and two orders of
magnitude cheaper than rebuilding the index, which is the comparison that
matters. (Nothing is `fsync`ed, so this is the cost of handing the bytes to
the OS rather than of a durable write.)

**hnswlib is about twice as quick at both**, and its file is ~8% smaller
(197 MB against 213 MB on Fashion-MNIST): it writes a fixed-size record per
element straight out of its own arena, where this library pickles, and it
stores no edge weights. Ours keeps them whenever `select` does — under
`select='distances'`, which stores none, the sizes match. Both are far below
the cost of building, so the difference is unlikely to be what you optimise.

**Only load files you trust.** `save`/`load` are built on `pickle`, and
unpickling any file runs code that the file chooses. That is pickle's design,
not something this library adds or can screen for: `load` sees the constructor
arguments only after the interpreter has already executed whatever the stream
asked for. Treat an index file exactly as you would a `.py` you were about to
run — the format is not a safe interchange format for data from elsewhere.

`load` is strict in both directions: it refuses a `metric=` for a file that
already has one (that would replace a metric known to be right with one it
cannot check), and refuses a metric of the wrong *kind* — the three forms store
their elements differently, so a callback cannot be attached to a graph built
for a built-in name. What it cannot check is that you passed the *same
distance*; supply a different one and the graph answers for a metric it was not
built for. That much is yours to get right.

Under the hood this is a plain pickle using
[`persistent_id`](https://docs.python.org/3/library/pickle.html#persistence-of-external-objects),
the stdlib's own mechanism for an object the loader must re-supply — there is
no bespoke format here to version. The pickle carries the vectors (or the
elements), the adjacency, the construction parameters and any soft deletions,
so it is self-contained but not small: expect roughly the in-memory size of the
index. There is no separate index-only/mmap format.

### The two classes

They take the same constructor arguments and answer the same calls, and both
accept all three metric forms. What only `HNSWFast` has is parallelism
(`update` across cores, `search_batch` across queries), numpy arrays out of the
batch API, and the `float32` vector store.

`HNSWPure` accepts `num_threads` and inserts serially, and rejects
`float32=True` having no typed store to narrow. It accepts
`keep_workspace=True` and ignores it, rather than rejecting it: `float32`
changes the numbers, so ignoring a request for it would mislead about results,
while `keep_workspace` only decides when memory is allocated and a class with
no threads runs the same program either way. It does call C metrics, through
`ctypes`, marshalling both operands on every distance. It is there so the
reference implementation can run the metric you
are actually indexing with — but that marshalling is a fixed per-call cost, so
whether it beats writing the distance in Python depends on how much work the
distance itself does. Against the plain Python loop a `jit_metric` wraps, the
pointer measured ~2.8x *slower* at `dim=8` and ~5x faster at `dim=64`. If you
just want a working index at low dimension, that original two-argument function
is still reachable on `.py_func_2arg`.

The two produce equivalent results only *distributionally* (ties are broken
differently), so compare them via k-NN recall rather than exact neighbor
lists.

## Performance

On the dense-vector case this is in the same class as the established C++
libraries — which is the point, since that case is the one they optimize.

**What the generality is and is not.** It is not that nobody else indexes
non-vector data: nmslib ships `jaccard_sparse`, `leven` and others, and is a
real competitor on anything in its catalogue. The difference is that a
catalogue is fixed at compile time and a Python callable is not. If your
dissimilarity is in nmslib's list, use nmslib and expect it to win; if it is a
business rule, a weighted ensemble, or anything you would otherwise have to
write C++ for, this library will run it — and, with `jit_metric`, run it at a
speed the section below quantifies.

### Against the other libraries

50k Gaussian vectors in 32 dimensions, m=16, efConstruction=200, k=10, 12
threads. Each figure is a ratio to this library, formed *within* a repetition
and shown with its range over three of them. A ratio whose range spans 1.0 is
not a difference; the resolution limit here is 1.06x on build and 1.07x on
query, measured by running the identical build twice as a control arm.

Both columns are **times**, relative to this library, so in both of them lower
is faster and 1.00x is us:

| library | build time | query time | recall@ef160 |
|---|---|---|--:|
| **this library** | 1.00x | 1.00x | .976 |
| hnswlib 0.8.0 | 1.03x [0.93-1.10] | 0.95x [0.94-1.06] | .980 |
| faiss `IndexHNSWFlat` | 1.06x [1.00-1.07] | 0.89x [0.88-0.93] | .979 |
| nmslib 2.1.2 (AVX2) | 1.44x [1.42-1.47] | 1.09x [1.06-1.12] | .980 |
| usearch | 1.36x [1.26-1.41] | 1.54x [1.49-1.54] | .978 |
| annoy (trees, not HNSW) | 0.17x [0.16-0.18] | 100x | 1.000 |

**hnswlib is indistinguishable from this library on both**, which is the
result that matters: the reference implementation of the algorithm, and we are
inside the noise of it. faiss queries ~11% faster. We build ~1.4x faster than
both nmslib and usearch, and query ~1.5x faster than usearch; nmslib's query
time is within a whisker of ours. annoy is a tree index included for
scale, and is the only arm whose parameters were not swept the way the others'
`ef` was.

Recall is identical to three decimals across every HNSW arm at matched `ef`,
which is the check that the comparison is like-for-like: the same graph
quality, differing only in how fast it is reached.

Two caveats. nmslib's PyPI wheel is built without SIMD, and warns so on
import; the row above is from a source build (`pip install --no-binary :all:
--no-build-isolation nmslib` — plain `--no-binary` fails on Python 3.14). And
the larger datasets we measured, at 100k x 128 and Fashion-MNIST, had control
spreads of 2.0x to 5.6x on this machine, so they are omitted rather than
reported: nothing there was resolvable.

### An arbitrary metric, measured

`kosarak` from ann-benchmarks: 74,962 sets of ints, mean size 55.6, Jaccard,
recall@10 against the published ground truth, M=16 efC=100. The comparison is
nmslib's `jaccard_sparse`, built from source so it has AVX2 — its purpose-built
C++ space against a Python function of ours.

Three ways to give this library the same distance. Every ratio below is
formed *inside* one run, from the minimum of seven alternating repetitions,
with a control arm — the same index timed a second time — so the resolution is
measured rather than hoped for. The two columns of ratios span ef=100 and
ef=400, that is recall 0.90 and 0.95:

| how | against `@jit_metric` | against nmslib | resolution |
|---|--:|--:|--:|
| plain Python callable | 7.5–8.1x slower | — | 1.02x |
| `@jit_metric`, 1 thread (pinned) | — | 1.18–1.20x slower | 1.002x |
| `@jit_metric`, 12 threads | — | 1.24–1.26x slower | 1.02–1.04x |

So **compiling the callable is worth about 8x**, and once compiled the gap to
a purpose-built C++ space is about **1.2x**, whether serial or threaded. That
is what writing your distance in Python instead of C++ costs.

Each arm queries through whichever API is faster for it: `search_batch` for
the compiled metric, a `search` loop for the Python callable, which serially
on this GIL build beat batching by ~6%. Matching them would widen the 8x
slightly, not narrow it.

For scale rather than for comparison, the serial pinned run measured 0.61
ms/query for `@jit_metric` against nmslib's 0.50 at recall 0.90, and 0.095
against 0.075 at twelve threads. Absolute figures move between runs on this
machine — the same compiled index measured 0.61 in one run and 0.73 in
another, after a long build had heated the box — which is exactly why the
ratios above are formed within a run and the absolutes are not compared
across them.

The serial rows are pinned to one core with `taskset`: this machine's cores
are not all the same speed, so an unpinned "1 thread" figure silently depends
on which one it landed on — measured 2.6x apart for the same code. See
[About these numbers](#about-these-numbers).

Variable-length objects reach a C metric through a fixed `double*` signature by
storing `[offset, length]` in the row and closing over one flat array, so the
vector store is 1.2 MB rather than the ~1.5 GB a padded matrix would need, and
the distance stays exact:

```python
FLAT = np.concatenate(sorted_sets)      # every set, back to back

@jit_metric
def jaccard(a, b):                      # each row is [offset, length]
    ao, an = int(a[0]), int(a[1])
    bo, bn = int(b[0]), int(b[1])
    i = j = inter = 0
    while i < an and j < bn:            # sorted-merge intersection
        x, y = FLAT[ao + i], FLAT[bo + j]
        if x == y:   inter += 1; i += 1; j += 1
        elif x < y:  i += 1
        else:        j += 1
    return 1.0 - inter / (an + bn - inter)
```

**When an index is worth building at all.** On an expensive metric the answer
depends on n, because brute force is linear and the graph is not. Jaccard over
random sets, ef=400, serial:

| n | brute force | this library | speedup |
|--:|--:|--:|--:|
| 5,000 | 5.04 ms/q | 4.68 ms/q | 1.1x |
| 20,000 | 22.25 ms/q | 9.50 ms/q | 2.3x |
| 80,000 | 87.72 ms/q | 12.81 ms/q | 6.8x |

Below ~10k elements an index barely pays for itself; by ~100k it is most of
the cost. (Recall drifts 0.955 -> 0.825 across those rows at fixed `ef`, so
part of the last speedup is bought with quality.)

### About these numbers

Every figure in this README comes from one throttling laptop, so the method
matters more than usual. What was done, and what it means for how much weight
a number can carry:

- **A control arm.** Every comparison includes the same build measured twice.
  The spread between those two runs is the resolution limit, and no effect
  smaller than it is claimed. This is not a formality: it caught a 10x
  "difference" that was two benchmark processes accidentally running at once,
  and a 3.8x one that was a metric interpolating across the steep part of a
  recall curve.
- **Interleaved, paired repetitions.** Arms alternate within a repetition and
  ratios are formed inside it, so drift in the machine is common-mode and
  cancels. Summarising each arm separately and dividing would not.
- **One process per measurement**, with the thread count pinned before any
  import and the observed value recorded, plus load average, CPU-seconds and
  timestamps on every row.
- **Pinning, where it decides anything.** The cores here are not all the same
  speed, so an unpinned single-threaded figure depends on which one it landed
  on — measured 2.6x apart for identical code. Serial comparisons use
  `taskset`; thread-count speedups are reported as measured, without being
  turned into an efficiency against a core count that would not mean much.

What this buys is roughly 3-13% resolution single-threaded and 6-7% on the
50k x 32 dataset at twelve threads. On the larger datasets at twelve threads it
buys nothing: identical builds differed by 2.0x to 5.6x, so those cells are
reported as unresolvable rather than given numbers. Redoing all of this on a
quiet, homogeneous machine would tighten every bound here, and is on the list.

**Reproducing any of this.** Everything above is produced by
[`benchmarks/`](benchmarks/) against the `anydist` you have installed —
`python run.py` then `python analyze.py`, with every comparison library
optional. The harness prints your machine's resolution before it prints any
comparison, so you can see which of its own numbers to believe.

## Author

Matteo Dell'Amico — `della@linux.it`

## License

BSD 3-clause; see the LICENSE file.
