Add IVF add and remove without the original data

.github/workflows/test.yaml

@@ -206,6 +206,7 @@ jobs:
           simd_filter+=':RaBitQKMeansTest.FinalModesKeepFlatLloydAssignmentsAcrossClusterCounts'
           simd_filter+=':IvfPersistenceTest.ReloadsRawAndQuantizedStorage'
           simd_filter+=':HnswSaveTest.SuccessfulSavePreservesSearchAfterLoading'
+          simd_filter+=':IvfAddTest.*:IvfRemoveTest.*'
           for cpu_model in hsw icl; do
             RABITQ_TEST_CPU="$cpu_model" "$RUNNER_TEMP/sde/sde64" "-$cpu_model" -- \
               build-install/tests/rabitq_tests --gtest_filter="$simd_filter"

DEVELOPMENT.md

@@ -354,6 +354,18 @@ compatibility test fixture or an explicit rejection path. IVF, HNSW, and Symphon
 formats and must each be reviewed. SymphonyQG includes a versioned quantized format plus a legacy
 raw-format fallback; preserve both unless a breaking change is explicitly requested.
 
+IVF `remove` stores a removed point as `f_add = +inf` inside the existing batch data, so it needs no
+discriminator. No index written before `remove` existed has an infinite `f_add`, so nothing in an old
+file is reinterpreted, and a file that has removals loads in release 0.5.0 too, where those points
+never appear in results. Keep this convention when touching IVF code:
+
+- Every estimator and scan path must turn `f_add = +inf` into an infinite estimated distance and an
+  infinite lower bound, and never into NaN. A NaN distance would enter the result buffer.
+- `fastscan::unpack_codes` must stay the exact inverse of `pack_codes`; `IVF::add` uses it to refill
+  the partial last batch of a cluster.
+- `IvfRemoveTest.ExcludesRemovedPointsOnEveryScanPath`, `IvfRemoveTest.SurvivesSaveLoadAndLaterAdds`,
+  and `FastScanPackingTest.UnpackInvertsPackIncludingTailsAndDoesNotOverrun` guard these rules.
+
 ### Change Python bindings
 
 Shared NumPy and string conversion helpers live in `python_bindings/bindings_common.hpp`. Register

docs/docs/index/ivf.md

@@ -156,3 +156,85 @@ remaining quantized bits or the stored raw vectors. Raw reranking computes
 squared L2 or `1 - dot(query, vector)` in the original coordinates; the search
 API is unchanged. Cluster selection and filtering remain approximate in both
 modes. Search returns the top `k` results after scanning the selected clusters.
+
+## Updating an Index
+A constructed or loaded index can gain and lose points without the original
+data. The index file format does not change.
+
+```c++
+void IVF::add(
+    const float* data,
+    size_t n,
+    const PID* cluster_ids = nullptr,
+    bool faster = false,
+    size_t num_threads = std::numeric_limits<size_t>::max()
+);
+
+size_t IVF::remove(const PID* ids_to_remove, size_t n);
+```
+
+- **data**, **n**: `n` new vectors, quantized like `construct` does. Vector `i`
+  receives the PID `max_elements() + i`, so PIDs stay dense.
+- **cluster_ids**: The cluster of each new vector, in `[0, cluster_num)`. When
+  it is `nullptr`, each vector goes to its nearest centroid, chosen the same way
+  a query is routed.
+- **faster**, **num_threads**: Same as in `construct`.
+- **ids_to_remove**: PIDs to remove. Every PID must be below `max_elements()`; nothing is
+  removed if one is not. `remove` returns how many were newly removed and can be
+  repeated safely.
+
+In Python:
+
+```python
+new_ids = index.add(vectors)                       # route to the nearest centroids
+new_ids = index.add(vectors, cluster_ids=labels)   # or choose the clusters
+removed = index.remove(new_ids[:10])
+```
+
+In C++, `construct` expects `max_elements()` rows once points have been added. It is
+not safe to call `add` or `remove` while another thread searches the same index.
+
+### Cost of `add`
+
+`add` builds the grown index next to the old one and swaps it in, so the index is
+unchanged if it throws. Every call copies the whole storage, so:
+
+- It needs memory for both copies while it runs.
+- Each call has a fixed cost that grows with the size of the whole index, on top of a
+  cost per new vector. **Add many vectors per call rather than one at a time.**
+
+On an index of 1,000,000 vectors with 128 dimensions, 1,024 clusters and 1 bit
+(8 threads, portable build, Intel i7-12700KF), adding 1,000 vectors in one call took
+about 5 ms, and adding the same vectors one call each took about 3 s (2.8 s and 3.7 s
+in two runs). For that 33 MB file, a process that had just loaded it held about
+65 MB, and peaked at about 100 MB during one `add` of 1,000 vectors.
+
+There is no spare capacity: every `add` moves the data.
+
+### Recall after `add`
+
+The centroids and the rotation are fixed when the index is constructed, so
+`add` never retrains them. If the added vectors come from a different distribution,
+or the index grows many times beyond its original size, recall at a given
+`nprobe` can drop. Rebuild from the original data with new centroids when
+recall matters more than the cost of a rebuild.
+
+### How removal is stored
+
+`remove` hides points from search but keeps their storage. Removed points still
+count in `max_elements()` and cannot be restored.
+
+A removed point has its `f_add` value set to `+inf`. That value is the constant
+term of the distance estimate, so the estimated distance and the lower bound of the
+point are `+inf` for every bit width, metric, and SIMD backend. The point never
+enters the result buffer and never triggers reranking. Any code that estimates
+IVF distances must keep this behavior and must never turn `+inf` into NaN.
+
+Search may return fewer than `k` results once points are removed: the missing slots
+hold `kPidMax` (`2**32 - 1` in Python) with an infinite distance, as when the probed
+clusters hold fewer than `k` points.
+
+The value lives in the ordinary batch data, so the file format does not change and
+removal survives `save` and `load`. No file written before `remove` existed has an
+infinite `f_add`, so nothing in an old file is reinterpreted. A file that has
+removals also loads in release 0.5.0, where the removed points never appear in results.

docs/docs/quick_start.md

@@ -97,7 +97,16 @@ ids, distances = index.search(queries, k=10, nprobe=5, high_accuracy=False)
 # Omit high_accuracy, or pass None, to use automatic selection.
 ```
 
-See the [IVF guide](index/ivf.md) for storage costs and persistence compatibility.
+An IVF index can also grow and shrink after it is built, without the original data:
+
+```python
+new_vectors = rng.standard_normal((50, 64)).astype(np.float32)
+new_ids = index.add(new_vectors)  # nearest-centroid routing; ids 500..549
+index.remove(new_ids[:10])        # hide ten of them from later searches
+```
+
+See the [IVF guide](index/ivf.md) for storage costs, persistence compatibility, and
+the limits of updating an index.
 
 The `metric` argument accepts `"l2"` and `"ip"` (also spelled
 `"innerproduct"`). To search by cosine similarity, normalize database and

include/rabitqlib/fastscan/fastscan.hpp

@@ -117,6 +117,43 @@ inline void pack_codes(
     }
 }
 
+/**
+ * @brief Inverse of pack_codes: recover the compact quantization codes of the first `num`
+ * vectors from their packed blocks. Padding lanes of a partial final batch are ignored.
+ *
+ * @param padded_dim dimension of quantized data (i.e., quantization code)
+ * @param blocks packed quantization code, as produced by pack_codes
+ * @param num   number of quantization code to recover
+ * @param quantization_code recovered codes, num * padded_dim / 8 bytes
+ */
+inline void unpack_codes(
+    size_t padded_dim, const uint8_t* blocks, size_t num, uint8_t* quantization_code
+) {
+    const size_t cols = padded_dim / 8;
+
+    for (size_t row = 0; row < num; row += kBatchSize) {
+        for (size_t i = 0; i < cols; ++i) {
+            for (size_t j = 0; j < 16; ++j) {
+                // blocks[j] holds the high nibbles of vector kPerm0[j] and the low nibbles
+                // of vector kPerm0[j] + 16 packed together; blocks[j + 16] holds the rest.
+                const size_t first = row + kPerm0[j];
+                const auto a = static_cast<uint8_t>(
+                    ((blocks[j] & 0x0F) << 4) | (blocks[j + 16] & 0x0F)
+                );
+                const auto b =
+                    static_cast<uint8_t>((blocks[j] & 0xF0) | (blocks[j + 16] >> 4));
+                if (first < num) {
+                    quantization_code[first * cols + i] = a;
+                }
+                if (first + 16 < num) {
+                    quantization_code[(first + 16) * cols + i] = b;
+                }
+            }
+            blocks += 32;
+        }
+    }
+}
+
 // Accumulate one batch into int32_t results. Throws if a result exceeds int32_t.
 // dim must be a positive multiple of 16.
 // Requires AVX2/FMA, AVX-512, or ARM NEON; otherwise throws std::runtime_error.