diff --git a/README.md b/README.md index 29e1351..75bef12 100644 --- a/README.md +++ b/README.md @@ -78,6 +78,12 @@ SAX summary for each insert using the breakpoints established during the initial routes each new summary into the live iSAX tree and splits full leaves without rebuilding the index. Its first insert copies a borrowed initial in-memory database into owned growable storage. +SAX encoding uses the ordered breakpoints for both tail intervals, including shifted and repeated +adaptive cuts. `setNormalized(false)` selects adaptive breakpoints without changing input values. +Python MESSI queries accept strided NumPy arrays (including `.fvecs` views); the binding copies +queries to contiguous float32 storage when needed and checks the index dimension before searching. +MESSI range refinement evaluates the full distance when a partial sum equals the inclusive radius. + `FaissFlat` is the streaming baseline: `faiss::IndexFlatL2::add()` appends into flat storage with no training step and no rebuild, so nothing is frozen at build time that could drift from the live data and the index stays exact after every insert. It supports `L2_SQUARED` only. diff --git a/docs/test-framework-guide.md b/docs/test-framework-guide.md index fb82926..dbe488a 100644 --- a/docs/test-framework-guide.md +++ b/docs/test-framework-guide.md @@ -28,3 +28,43 @@ Follow the parameterized test pattern used throughout the test suite for BruteFo - Tests should be parameterized to run with different configurations (k-values, thread counts, datasets) - Include proper ground truth validation - Follow the existing naming conventions and code structure + +## MESSI exact top-k regressions + +`test_Messi_TopK_Distinct` compares every returned neighbor against exhaustive +brute-force ground truth, checks distinct IDs, and recomputes the squared-L2 +distance of each returned vector. Missing neighbors are failures; a recall +threshold below 100% is not sufficient for exact search. Integer-distance tie +fixtures allow any distinct IDs with the exact expected distances. + +The parameterized cases cover Gaussian and adaptive breakpoints, one and four +query workers, k=1/10/100/1000, repeated queries, and streaming updates. The +inclusive-range regression checks zero and nonzero radii, exact boundary hits and +partial sums equal to the radius before later nonzero chunks. The adaptive MESSI +case in `test_equidepth` checks all k results with a small floating-point +tolerance rather than only checking finite output. Simple CI +runs these regressions through its complete CTest suite. + +`test_sax_breakpoints` covers positive/negative shifted cuts, repeated cuts, exact +boundaries and both tail bins. It also checks scalar and SIMD record lower bounds +against exact squared L2. The tests fail with the old sign-based SAX fallback. + +With `BUILD_PYTHON=ON`, `MessiPythonLayouts` additionally checks exact top-k and +inclusive range results for contiguous, `.fvecs`-style, row-strided, Fortran and +negative-stride queries, before and after streaming inserts. NumPy must be +installed for this Python test. Python builds with Coconut disabled omit its binding. + +For the original unnormalized SIFT1M data, build Python bindings and install +`numpy` and `faiss-cpu`, then run the optional manual audit: + +```bash +PYTHONPATH=. python tests/check_messi_sift.py --dataset /path/to/sift --output sift-report.json +``` + +The audit uses 100 queries, k=1/10/100/1000, Gaussian/adaptive modes, one/four +search workers, repeated queries, static 20K/1M builds and streaming 750K→875K→1M. +It checks every distance against exact FAISS IndexFlatL2, recomputes returned +ID/distance pairs, rejects duplicates and validates the full-dataset reference +against the published SIFT ground truth. Distances of raw integer descriptors +are compared exactly; equally distant alternative IDs are allowed. The dataset +and FAISS are not required for the regular C++ tests. diff --git a/lib/algos/Messi.cpp b/lib/algos/Messi.cpp index 14e7cd3..67fb50c 100644 --- a/lib/algos/Messi.cpp +++ b/lib/algos/Messi.cpp @@ -128,7 +128,7 @@ namespace daisy float dist = ts_euclidean_distance_SIMD(query, node->buffer->full_ts_buffer[i], index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]); - if (dist <= pq_bsf->knn[pq_bsf->k - 1]) + if (dist < pq_bsf->knn[pq_bsf->k - 1]) { file_position_type pos = 0; if (node->buffer->full_position_buffer != nullptr && @@ -137,7 +137,7 @@ namespace daisy pos = *node->buffer->full_position_buffer[i] / index->settings->timeseries_size; } pthread_rwlock_wrlock(lock_queue); - pqueue_bsf_insert(pq_bsf, dist, static_cast(pos), node); + pqueue_bsf_insert_no_dup_position(pq_bsf, dist, static_cast(pos), node); pthread_rwlock_unlock(lock_queue); } } @@ -147,7 +147,7 @@ namespace daisy float dist = ts_euclidean_distance_SIMD(query, node->buffer->tmp_full_ts_buffer[i], index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]); - if (dist <= pq_bsf->knn[pq_bsf->k - 1]) + if (dist < pq_bsf->knn[pq_bsf->k - 1]) { file_position_type pos = 0; if (node->buffer->tmp_full_position_buffer != nullptr && @@ -156,7 +156,7 @@ namespace daisy pos = *node->buffer->tmp_full_position_buffer[i] / index->settings->timeseries_size; } pthread_rwlock_wrlock(lock_queue); - pqueue_bsf_insert(pq_bsf, dist, static_cast(pos), node); + pqueue_bsf_insert_no_dup_position(pq_bsf, dist, static_cast(pos), node); pthread_rwlock_unlock(lock_queue); } } @@ -174,11 +174,11 @@ namespace daisy float dist = ts_euclidean_distance_SIMD(query, &(rawfile[*node->buffer->partial_position_buffer[i]]), index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]); - if (dist <= pq_bsf->knn[pq_bsf->k - 1]) + if (dist < pq_bsf->knn[pq_bsf->k - 1]) { pthread_rwlock_wrlock(lock_queue); - pqueue_bsf_insert(pq_bsf, dist, *node->buffer->partial_position_buffer[i] / index->settings->timeseries_size, node); + pqueue_bsf_insert_no_dup_position(pq_bsf, dist, *node->buffer->partial_position_buffer[i] / index->settings->timeseries_size, node); pthread_rwlock_unlock(lock_queue); } @@ -229,7 +229,7 @@ namespace daisy if (dist <= pq_bsf->knn[pq_bsf->k - 1]) { pthread_rwlock_wrlock(lock_queue); - pqueue_bsf_insert(pq_bsf, dist, *node->buffer->partial_position_buffer[i] / index->settings->timeseries_size, node); + pqueue_bsf_insert_no_dup_position(pq_bsf, dist, *node->buffer->partial_position_buffer[i] / index->settings->timeseries_size, node); pthread_rwlock_unlock(lock_queue); } } @@ -496,10 +496,13 @@ namespace daisy if (node == NULL || node->buffer == NULL) return; + // The bounded SIMD kernel stops at >= bound. Inclusive ranges must + // continue past partial sums equal to r before accepting a full distance. + const float evaluation_bound = std::nextafter(r, std::numeric_limits::infinity()); for (int i = 0; i < node->buffer->full_buffer_size; i++) { float dist = ts_euclidean_distance_SIMD(query, node->buffer->full_ts_buffer[i], - index->settings->timeseries_size, r); + index->settings->timeseries_size, evaluation_bound); if (dist <= r) { file_position_type pos = 0; @@ -515,7 +518,7 @@ namespace daisy for (int i = 0; i < node->buffer->tmp_full_buffer_size; i++) { float dist = ts_euclidean_distance_SIMD(query, node->buffer->tmp_full_ts_buffer[i], - index->settings->timeseries_size, r); + index->settings->timeseries_size, evaluation_bound); if (dist <= r) { file_position_type pos = 0; @@ -540,7 +543,7 @@ namespace daisy if (distmin <= r) { float dist = ts_euclidean_distance_SIMD(query, &(rawfile[*node->buffer->partial_position_buffer[i]]), - index->settings->timeseries_size, r); + index->settings->timeseries_size, evaluation_bound); if (dist <= r) { pthread_rwlock_wrlock(lock); @@ -558,6 +561,8 @@ namespace daisy ts_type *paa = wd->paa; ts_type *ts = wd->ts; float r = wd->r; + // The shared top-k enqueue helper uses <; range pruning is inclusive. + const float pruning_bound = std::nextafter(r, std::numeric_limits::infinity()); int n_pqueue = wd->n_pqueue; float *rawfile = wd->rawfile; int tnumber = rand() % n_pqueue; @@ -570,7 +575,7 @@ namespace daisy int current_root_node_number = __sync_fetch_and_add(wd->node_counter, 1); if (current_root_node_number >= wd->amountnode) break; - insert_tree_node_m_hybridpqueue(paa, wd->nodelist[current_root_node_number], index, r, + insert_tree_node_m_hybridpqueue(paa, wd->nodelist[current_root_node_number], index, pruning_bound, wd->allpq, wd->alllock, &tnumber, n_pqueue); } diff --git a/lib/isax/SAX.cpp b/lib/isax/SAX.cpp index bc537bc..be8fe16 100644 --- a/lib/isax/SAX.cpp +++ b/lib/isax/SAX.cpp @@ -1,5 +1,6 @@ #include #include +#include #include "SAX.hpp" #include "../utils/SIMD.hpp" @@ -86,22 +87,13 @@ namespace daisy // TO (c - 1) * (c - 2) / 2 + c - 1 int offset = ((cardinality - 1) * (cardinality - 2)) / 2; - int si; - for (si = 0; si < segments; si++) - { - sax_out[si] = 0; - - // First object = daisy_active_breakpoints[offset] - // Last object = daisy_active_breakpoints[offset + cardinality - 2] - // Size of sub-array = cardinality - 1 - - float *res = (float *)bsearch(&paa[si], &daisy_active_breakpoints[offset], cardinality - 1, - sizeof(ts_type), compare); - if (res != NULL) - sax_out[si] = (int)(res - &daisy_active_breakpoints[offset]); - else if (paa[si] > 0) - sax_out[si] = cardinality - 1; - } + // Each interval is (previous_cut, cut], with unbounded tail bins. + // Searching the cuts directly also handles shifted/repeated adaptive + // cuts; neither the preceding cardinality nor the sign is an endpoint. + const float *cuts = daisy_active_breakpoints + offset; + for (int si = 0; si < segments; si++) + sax_out[si] = static_cast( + std::lower_bound(cuts, cuts + cardinality - 1, paa[si]) - cuts); free(paa); return SUCCESS; @@ -126,22 +118,10 @@ namespace daisy { int offset = ((cardinality - 1) * (cardinality - 2)) / 2; - int si; - for (si = 0; si < segments; si++) - { - sax[si] = 0; - - // First object = daisy_active_breakpoints[offset] - // Last object = daisy_active_breakpoints[offset + cardinality - 2] - // Size of sub-array = cardinality - 1 - - float *res = (float *)bsearch(&paa[si], &daisy_active_breakpoints[offset], cardinality - 1, - sizeof(ts_type), compare); - if (res != NULL) - sax[si] = (int)(res - &daisy_active_breakpoints[offset]); - else if (paa[si] > 0) - sax[si] = cardinality - 1; - } + const float *cuts = daisy_active_breakpoints + offset; + for (int si = 0; si < segments; si++) + sax[si] = static_cast( + std::lower_bound(cuts, cuts + cardinality - 1, paa[si]) - cuts); return SUCCESS; } diff --git a/lib/isax/iSAXPqueue.cpp b/lib/isax/iSAXPqueue.cpp index 71c0ce9..545a36c 100644 --- a/lib/isax/iSAXPqueue.cpp +++ b/lib/isax/iSAXPqueue.cpp @@ -498,11 +498,26 @@ namespace daisy { int i, j; - // Fast duplicate check: if position already present, skip + // If the position is already queued, keep only its best distance. An + // earlier entry may carry an early-abandoned (over-estimated) value, so + // a better distance replaces it instead of being dropped. for (i = 0; i < q->k; i++) { if (q->position[i] == position) - return; + { + if (data >= q->knn[i]) + return; + for (j = i; j < q->k - 1; j++) + { + q->knn[j] = q->knn[j + 1]; + q->position[j] = q->position[j + 1]; + q->node[j] = q->node[j + 1]; + } + q->knn[q->k - 1] = FLT_MAX; + q->position[q->k - 1] = -1; + q->node[q->k - 1] = nullptr; + break; + } } for (i = 0; i < q->k; i++) diff --git a/lib/isax/iSAXSearch.cpp b/lib/isax/iSAXSearch.cpp index 3090a96..23c783c 100644 --- a/lib/isax/iSAXSearch.cpp +++ b/lib/isax/iSAXSearch.cpp @@ -22,7 +22,7 @@ namespace daisy { float dist = ts_euclidean_distance_SIMD(query, node->buffer->full_ts_buffer[i], index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]); - if (dist <= pq_bsf->knn[pq_bsf->k - 1]) + if (dist < pq_bsf->knn[pq_bsf->k - 1]) { file_position_type pos = 0; if (node->buffer->full_position_buffer != nullptr && @@ -36,7 +36,7 @@ namespace daisy { float dist = ts_euclidean_distance_SIMD(query, node->buffer->tmp_full_ts_buffer[i], index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]); - if (dist <= pq_bsf->knn[pq_bsf->k - 1]) + if (dist < pq_bsf->knn[pq_bsf->k - 1]) { file_position_type pos = 0; if (node->buffer->tmp_full_position_buffer != nullptr && @@ -53,7 +53,7 @@ namespace daisy float dist = ts_euclidean_distance_SIMD(query, &(rawfile[*node->buffer->partial_position_buffer[i]]), index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]); - if (dist <= pq_bsf->knn[pq_bsf->k - 1]) + if (dist < pq_bsf->knn[pq_bsf->k - 1]) { file_position_type pos = *node->buffer->partial_position_buffer[i] / index->settings->timeseries_size; pqueue_bsf_insert_no_dup_position(pq_bsf, dist, pos, node); diff --git a/pybinds/setup.cpp b/pybinds/setup.cpp index ca227a6..2c1636c 100644 --- a/pybinds/setup.cpp +++ b/pybinds/setup.cpp @@ -12,7 +12,9 @@ #include "../lib/distance_computers/DistanceComputer.hpp" #include "../lib/algos/Bruteforce.hpp" #include "../lib/algos/LbBruteforce.hpp" +#ifdef COCONUT_ENABLED #include "../lib/algos/Coconut.hpp" +#endif #include "../lib/algos/Messi.hpp" #if ODYSSEY_MPI #include "../lib/algos/hodyssey/Odyssey.hpp" @@ -377,6 +379,7 @@ PYBIND11_MODULE(_core, m) return pybind11::make_tuple(I, D); }, "Search using SearchConfig (top-k or range) and return (indices, distances)"); ////// COCONUT ////// +#ifdef COCONUT_ENABLED pybind11::class_(m, "Coconut", "COCONUT sortable-SAX index (static bottom-up build + streaming insert)") .def(pybind11::init(), "Create a new Coconut with the given distance metric") @@ -447,6 +450,8 @@ PYBIND11_MODULE(_core, m) self.searchIndex(static_cast(query_buf.ptr), n_query, config, I, D); return pybind11::make_tuple(I, D); }, "Search using SearchConfig (top-k or range) and return (indices, distances)"); +#endif + ////// MESSI ////// pybind11::class_(m, "Messi", "MESSI (Multi-Queue Efficient SAX Similarity Index) algorithm for time series similarity search") // Constructor @@ -534,16 +539,22 @@ PYBIND11_MODULE(_core, m) "Incrementally insert a batch of series into the live MESSI index") // Search the index with query array and return top-k results - .def("searchIndex", [](daisy::Messi &self, pybind11::array_t query, daisy::idx_t k) + .def("searchIndex", [](daisy::Messi &self, + pybind11::array_t query, + daisy::idx_t k) { pybind11::buffer_info query_buf = query.request(); if (query_buf.ndim != 2) throw std::runtime_error("Query array must be 2D"); + if (self.getDim() == 0 || self.getNDatabase() == 0) + throw std::runtime_error("MESSI index must be built before searching"); + if (static_cast(query_buf.shape[1]) != self.getDim()) + throw std::runtime_error("Query dimension does not match the index dimension"); if (k <= 0) throw std::runtime_error("k must be positive"); const daisy::idx_t n_query = query_buf.shape[0]; - const daisy::idx_t dim = query_buf.shape[1]; std::vector indices(n_query * k); std::vector distances(n_query * k); @@ -554,11 +565,20 @@ PYBIND11_MODULE(_core, m) pybind11::array_t({n_query, k}, indices.data()), pybind11::array_t({n_query, k}, distances.data()) ); }, "Search the MESSI index using queries and return (indices, distances)") - .def("searchIndex", [](daisy::Messi &self, pybind11::array_t query, daisy::SearchConfig config) + .def("searchIndex", [](daisy::Messi &self, + pybind11::array_t query, + daisy::SearchConfig config) { pybind11::buffer_info query_buf = query.request(); if (query_buf.ndim != 2) throw std::runtime_error("Query array must be 2D"); + if (self.getDim() == 0 || self.getNDatabase() == 0) + throw std::runtime_error("MESSI index must be built before searching"); + if (static_cast(query_buf.shape[1]) != self.getDim()) + throw std::runtime_error("Query dimension does not match the index dimension"); + if (config.type == daisy::QueryType::TOP_K && config.k == 0) + throw std::runtime_error("k must be positive"); daisy::idx_t n_query = query_buf.shape[0]; std::vector> I; std::vector> D; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8dad2f7..7708c70 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -370,6 +370,29 @@ gtest_discover_tests(test_LbBruteforce_Streaming WORKING_DIRECTORY ${CMAKE_SOURC message(STATUS "Tests discovered for test_Messi_L2Square.") endif() +# ////// MESSI top-k distinct results ////// +add_executable( + test_Messi_TopK_Distinct + test_Messi_TopK_Distinct.cpp +) +target_link_libraries( + test_Messi_TopK_Distinct + PRIVATE + GTest::gtest_main + dino_lib + commons_lib + stdc++fs +) +target_include_directories(test_Messi_TopK_Distinct + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons +) +gtest_discover_tests( + test_Messi_TopK_Distinct + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} +) + # ////// MESSI streaming ////// add_executable( test_Messi_Streaming @@ -1114,3 +1137,18 @@ gtest_discover_tests( test_equidepth WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} ) + +# SAX encoding must keep lower bounds valid for shifted/repeated adaptive cuts. +add_executable(test_sax_breakpoints test_sax_breakpoints.cpp) +target_link_libraries(test_sax_breakpoints PRIVATE GTest::gtest_main dino_lib commons_lib stdc++fs) +target_include_directories(test_sax_breakpoints PRIVATE ${CMAKE_SOURCE_DIR}/lib) +gtest_discover_tests(test_sax_breakpoints WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}) + +if(BUILD_PYTHON) + find_package(Python3 COMPONENTS Interpreter REQUIRED) + add_test(NAME MessiPythonLayouts + COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_messi_query_layout.py) + set_tests_properties(MessiPythonLayouts PROPERTIES + ENVIRONMENT "PYTHONPATH=${CMAKE_SOURCE_DIR}" + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}) +endif() diff --git a/tests/check_messi_sift.py b/tests/check_messi_sift.py new file mode 100644 index 0000000..7bfd33b --- /dev/null +++ b/tests/check_messi_sift.py @@ -0,0 +1,102 @@ +"""Manual SIFT exactness audit. Requires numpy, faiss-cpu and built daisy bindings. + +python tests/check_messi_sift.py --dataset /path/to/sift --output report.json +Dataset must contain sift_base.fvecs, sift_query.fvecs, sift_groundtruth.ivecs. +Uses raw descriptors, exact squared L2, and permits arbitrary IDs at exact ties. +""" +import argparse +import json +from pathlib import Path +import sys +import numpy as np +import faiss +import daisy + + +def vectors(path, dtype): + raw = np.memmap(path, dtype=np.int32, mode="r") + width = int(raw[0]) + rows = raw.reshape(-1, width + 1) + if not np.all(rows[:, 0] == width): + raise ValueError(f"Invalid dimension headers: {path}") + return rows.view(dtype)[:, 1:] + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--dataset", type=Path, required=True) + parser.add_argument("--output", type=Path, required=True) + parser.add_argument("--queries", type=int, default=100) + args = parser.parse_args() + base = np.ascontiguousarray(vectors(args.dataset / "sift_base.fvecs", np.float32)) + query = vectors(args.dataset / "sift_query.fvecs", np.float32)[:args.queries] + published = vectors(args.dataset / "sift_groundtruth.ivecs", np.int32)[:args.queries] + if base.shape != (1000000, 128): + raise ValueError(f"Expected original SIFT1M, got {base.shape}") + # Use direct L2 rather than cancellation-prone norm/dot-product evaluation. + faiss.omp_set_num_threads(4) + faiss.cvar.distance_compute_blas_threshold = 1000000000 + report = {"dataset": str(args.dataset.resolve()), "normalized_input": False, + "query_shape": list(query.shape), "query_strides": list(query.strides), + "module": daisy._core.__file__, "cases": []} + reference = {} + + def save(): + args.output.parent.mkdir(parents=True, exist_ok=True) + args.output.write_text(json.dumps(report, indent=2) + "\n") + + def expected(size): + if size not in reference: + index = faiss.IndexFlatL2(base.shape[1]) + index.add(base[:size]) + distances, _ = index.search(np.ascontiguousarray(query), 1000) + reference[size] = distances + if size == len(base): + stored = ((base[published].astype(np.float64) - + query.astype(np.float64)[:, None, :]) ** 2).sum(axis=2) + np.testing.assert_array_equal(distances[:, :published.shape[1]], stored) + report["published_groundtruth_matches"] = True + print(f"reference ready: {size}", flush=True) + return reference[size] + + # Small static reproduction, full static index, and large streaming index. + for stages in ((20000,), (1000000,), (750000, 875000, 1000000)): + for normalized in (False, True): + for workers in (1, 4): + index = daisy.Messi(daisy.DistanceType.L2_SQUARED) + index.setNormalized(normalized) + index.setSearchWorkers(workers) + index.setIndexWorkers(2) + initial = np.ascontiguousarray(base[:stages[0]]) + index.buildIndex(initial) + previous = stages[0] + for size in stages: + if size != previous: + index.insertBatch(base[previous:size]) + previous = size + gt = expected(size) + for k in (1, 10, 100, 1000): + for repeat in range(2): + # Deliberately retain fvecs header strides at the Python API. + ids, distances = index.searchIndex(query, k) + ids = np.asarray(ids, dtype=np.int64) + physical = ((base[ids].astype(np.float64) - + query.astype(np.float64)[:, None, :]) ** 2).sum(axis=2) + case = {"initial": stages[0], "size": size, + "mode": "gaussian" if normalized else "adaptive", + "workers": workers, "k": k, "repeat": repeat, + "wrong_distances": int(np.count_nonzero(distances != gt[:, :k])), + "invalid_id_distance_pairs": int(np.count_nonzero(distances != physical)), + "duplicate_rows": sum(len(set(row)) != k for row in ids)} + report["cases"].append(case) + save() + print(json.dumps(case), flush=True) + if any(case[key] for key in ("wrong_distances", "invalid_id_distance_pairs", "duplicate_rows")): + return 1 + report["passed"] = True + save() + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/test_Messi_TopK_Distinct.cpp b/tests/test_Messi_TopK_Distinct.cpp new file mode 100644 index 0000000..fb1b156 --- /dev/null +++ b/tests/test_Messi_TopK_Distinct.cpp @@ -0,0 +1,281 @@ +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +// Regression test: MESSI's approximate phase scans one leaf to seed the +// best-so-far queue, then the parallel phase scans every leaf again. Without +// de-duplication by position the seeded leaf's series were inserted twice, so +// the k-NN list contained the same series several times and displaced real +// neighbours. +namespace +{ + constexpr daisy::idx_t kRows = 20000; + constexpr daisy::idx_t kDim = 128; + constexpr daisy::idx_t kQueries = 50; + constexpr daisy::idx_t kK = 10; + + // 50 Gaussian clusters with isotropic noise: dense neighbourhoods, so the + // leaf scanned by the approximate phase holds members of the final top-k. + std::vector clusteredData(daisy::idx_t rows, unsigned seed) + { + std::mt19937 centre_rng(42); + std::normal_distribution normal; + std::vector centres(50 * kDim); + for (auto &v : centres) + v = 6.0F * normal(centre_rng); + std::mt19937 rng(seed); + std::uniform_int_distribution pick(0, 49); + std::vector data(static_cast(rows) * kDim); + for (daisy::idx_t r = 0; r < rows; ++r) + { + const float *centre = centres.data() + pick(rng) * kDim; + for (daisy::idx_t d = 0; d < kDim; ++d) + data[static_cast(r) * kDim + d] = centre[d] + normal(rng); + } + return data; + } +} + +TEST(MessiTopKDistinctTest, ReturnsDistinctIdsAndMatchesBruteForce) +{ + auto database = clusteredData(kRows, 1); + auto queries = clusteredData(kQueries, 2); + + daisy::BruteForceSearch brute(daisy::DistanceType::L2_SQUARED); + brute.buildIndex(database.data(), kRows, kDim); + std::vector expected_ids(kQueries * kK); + std::vector expected_distances(kQueries * kK); + brute.searchIndex(queries.data(), kQueries, kK, expected_ids.data(), + expected_distances.data()); + + daisy::MessiConfig config; + config.search_workers = 4; + config.index_workers = 2; + daisy::Messi messi(daisy::DistanceType::L2_SQUARED, config); + messi.setNormalized(false); + messi.buildIndex(database.data(), kRows, kDim); + std::vector ids(kQueries * kK); + std::vector distances(kQueries * kK); + messi.searchIndex(queries.data(), kQueries, kK, ids.data(), distances.data()); + + for (daisy::idx_t q = 0; q < kQueries; ++q) + { + std::set returned(ids.begin() + q * kK, ids.begin() + (q + 1) * kK); + EXPECT_EQ(returned.size(), static_cast(kK)) + << "query " << q << " returned duplicate ids"; + EXPECT_EQ(ids[q * kK], expected_ids[q * kK]) << "nearest neighbour differs, query " << q; + for (daisy::idx_t j = 0; j < kK; ++j) + { + EXPECT_EQ(returned.count(expected_ids[q * kK + j]), 1u) + << "missing exact neighbor, query " << q << ", rank " << j; + ASSERT_LT(ids[q * kK + j], kRows); + double actual = 0.0; + for (daisy::idx_t d = 0; d < kDim; ++d) + { + double delta = static_cast(database[ids[q * kK + j] * kDim + d]) - queries[q * kDim + d]; + actual += delta * delta; + } + EXPECT_NEAR(distances[q * kK + j], actual, 1e-5 + 1e-5 * std::fabs(actual)) + << "returned ID/distance pair differs, query " << q << ", rank " << j; + EXPECT_NEAR(distances[q * kK + j], expected_distances[q * kK + j], + 1e-5 + 1e-5 * std::fabs(expected_distances[q * kK + j])) + << "exact top-k distance differs, query " << q << ", rank " << j; + } + } +} + +class MessiExactnessTest : public ::testing::TestWithParam> {}; + +TEST_P(MessiExactnessTest, RemainsExactAcrossStreamingStagesAndRepeatedQueries) +{ + constexpr daisy::idx_t dim = 64, rows = 4096, nquery = 12; + const auto [normalized, workers] = GetParam(); + std::mt19937 rng(731); + std::normal_distribution normal; + std::vector database(rows * dim), query(nquery * dim); + auto fill = [&](std::vector &values, daisy::idx_t n) { + for (daisy::idx_t r = 0; r < n; ++r) + { + double mean = 0.0, variance = 0.0; + for (daisy::idx_t d = 0; d < dim; ++d) + { + float value = normal(rng); + values[r * dim + d] = normalized ? value : 75.0f + 8.0f * value; + mean += values[r * dim + d]; + } + if (normalized) + { + mean /= dim; + for (daisy::idx_t d = 0; d < dim; ++d) + variance += std::pow(values[r * dim + d] - mean, 2); + const double scale = std::sqrt(variance / dim); + for (daisy::idx_t d = 0; d < dim; ++d) + values[r * dim + d] = (values[r * dim + d] - mean) / scale; + } + } + }; + fill(database, rows); + fill(query, nquery); + // Include exact matches in the initial tree and both appended batches. + for (daisy::idx_t q = 0; q < 3; ++q) + std::copy_n(database.data() + (q == 0 ? 17 : q == 1 ? 3200 : 4000) * dim, + dim, query.data() + q * dim); + daisy::MessiConfig config; + config.index_workers = 2; + config.search_workers = workers; + daisy::Messi messi(daisy::DistanceType::L2_SQUARED, config); + messi.setNormalized(normalized); + messi.buildIndex(database.data(), 3072, dim); + daisy::idx_t previous = 3072; + for (daisy::idx_t n : {3072ULL, 3584ULL, 4096ULL}) + { + if (n > previous) + messi.insertBatch(database.data() + previous * dim, n - previous); + previous = n; + daisy::BruteForceSearch brute(daisy::DistanceType::L2_SQUARED); + brute.buildIndex(database.data(), n, dim); + for (daisy::idx_t k : {1ULL, 10ULL, 100ULL, 1000ULL}) + { + std::vector expected_ids(nquery * k), ids(nquery * k); + std::vector expected_distances(nquery * k), distances(nquery * k); + brute.searchIndex(query.data(), nquery, k, expected_ids.data(), expected_distances.data()); + for (int repeat = 0; repeat < 2; ++repeat) + { + SCOPED_TRACE(::testing::Message() << "normalized=" << normalized << ", workers=" << workers + << ", n=" << n << ", k=" << k << ", repeat=" << repeat); + messi.searchIndex(query.data(), nquery, k, ids.data(), distances.data()); + for (daisy::idx_t q = 0; q < nquery; ++q) + { + std::set returned(ids.begin() + q * k, ids.begin() + (q + 1) * k); + EXPECT_EQ(returned.size(), k) << "duplicate IDs, query " << q; + for (daisy::idx_t j = 0; j < k; ++j) + { + EXPECT_EQ(returned.count(expected_ids[q * k + j]), 1u) + << "missing exact neighbor, query " << q << ", rank " << j; + ASSERT_LT(ids[q * k + j], n); + double actual = 0.0; + for (daisy::idx_t d = 0; d < dim; ++d) + { + double delta = static_cast(database[ids[q * k + j] * dim + d]) - query[q * dim + d]; + actual += delta * delta; + } + EXPECT_NEAR(distances[q * k + j], actual, 1e-5 + 1e-5 * actual) + << "returned ID/distance mismatch, query " << q << ", rank " << j; + EXPECT_NEAR(distances[q * k + j], expected_distances[q * k + j], + 1e-5 + 1e-5 * expected_distances[q * k + j]) + << "exact distance mismatch, query " << q << ", rank " << j; + } + } + } + } + } +} + +INSTANTIATE_TEST_SUITE_P(BreakpointsAndWorkers, MessiExactnessTest, + ::testing::Combine(::testing::Bool(), ::testing::Values(1, 4))); + +TEST(MessiTopKDistinctTest, ExactTiesAllowAnyDistinctEquallyDistantNeighbors) +{ + constexpr daisy::idx_t dim = 64, rows = 144, nquery = 4; + std::vector database(rows * dim), query(nquery * dim); + for (daisy::idx_t r = 0; r < rows; ++r) + std::fill_n(database.data() + r * dim, dim, r < 40 || r >= 128 ? 0.0f : r < 80 ? 1.0f : 2.0f); + for (daisy::idx_t q = 0; q < nquery; ++q) + std::fill_n(query.data() + q * dim, dim, static_cast(q)); + daisy::MessiConfig config; + config.index_workers = 2; + config.search_workers = 4; + daisy::Messi messi(daisy::DistanceType::L2_SQUARED, config); + messi.setNormalized(false); + messi.buildIndex(database.data(), 96, dim); + for (daisy::idx_t n : {96ULL, rows}) + { + if (n == rows) + messi.insertBatch(database.data() + 96 * dim, rows - 96); + daisy::BruteForceSearch brute(daisy::DistanceType::L2_SQUARED); + brute.buildIndex(database.data(), n, dim); + for (daisy::idx_t k : {1ULL, 10ULL, 80ULL}) + { + std::vector expected_ids(nquery * k), ids(nquery * k); + std::vector expected_distances(nquery * k), distances(nquery * k); + brute.searchIndex(query.data(), nquery, k, expected_ids.data(), expected_distances.data()); + messi.searchIndex(query.data(), nquery, k, ids.data(), distances.data()); + for (daisy::idx_t q = 0; q < nquery; ++q) + { + std::set returned(ids.begin() + q * k, ids.begin() + (q + 1) * k); + EXPECT_EQ(returned.size(), k); + for (daisy::idx_t j = 0; j < k; ++j) + { + ASSERT_LT(ids[q * k + j], n); + float actual = 0.0f; + for (daisy::idx_t d = 0; d < dim; ++d) + { + float delta = database[ids[q * k + j] * dim + d] - query[q * dim + d]; + actual += delta * delta; + } + EXPECT_FLOAT_EQ(distances[q * k + j], actual); + EXPECT_FLOAT_EQ(distances[q * k + j], expected_distances[q * k + j]); + } + } + } + } +} + +TEST_P(MessiExactnessTest, InclusiveRangeRejectsPartialSumAtRadius) +{ + constexpr daisy::idx_t dim = 128, rows = 32; + const auto [normalized, workers] = GetParam(); + std::vector database(rows * dim, 0.0f), query(dim, 0.0f); + // Exact boundary hit: first SIMD chunk contributes 8, all others zero. + std::fill_n(database.data() + dim, 8, 1.0f); + // False boundary hit if refinement abandons after the first chunk. + std::fill_n(database.data() + 2 * dim, 8, 1.0f); + std::fill_n(database.data() + 3 * dim - 8, 8, 1.0f); + // Radius zero must also inspect chunks following an all-zero prefix. + database[3 * dim + dim - 1] = 1.0f; + daisy::MessiConfig config; + config.index_workers = 1; + config.search_workers = workers; + daisy::Messi messi(daisy::DistanceType::L2_SQUARED, config); + messi.setNormalized(normalized); + messi.buildIndex(database.data(), 16, dim); + for (daisy::idx_t size : {16, 32}) { + if (size == rows) + messi.insertBatch(database.data() + 16 * dim, 16); + for (float radius : {0.0f, 8.0f, 16.0f}) { + daisy::SearchConfig search; + search.type = daisy::QueryType::RANGE; + search.r = radius; + std::vector> ids; + std::vector> distances; + messi.searchIndex(query.data(), 1, search, ids, distances); + std::set expected; + for (daisy::idx_t id = 0; id < size; ++id) { + float exact = 0; + for (daisy::idx_t d = 0; d < dim; ++d) + exact += database[id * dim + d] * database[id * dim + d]; + if (exact <= radius) expected.insert(id); + } + ASSERT_EQ(ids.size(), 1u); + ASSERT_EQ(distances.size(), 1u); + EXPECT_EQ(std::set(ids[0].begin(), ids[0].end()), expected); + ASSERT_EQ(ids[0].size(), expected.size()); + ASSERT_EQ(ids[0].size(), distances[0].size()); + for (size_t j = 0; j < ids[0].size(); ++j) { + ASSERT_LT(ids[0][j], size); + float exact = 0; + for (daisy::idx_t d = 0; d < dim; ++d) + exact += database[ids[0][j] * dim + d] * database[ids[0][j] * dim + d]; + EXPECT_EQ(distances[0][j], exact); + } + } + } +} diff --git a/tests/test_equidepth.cpp b/tests/test_equidepth.cpp index 3c11a40..0fbc0f2 100644 --- a/tests/test_equidepth.cpp +++ b/tests/test_equidepth.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include @@ -42,8 +43,8 @@ namespace return data; } - constexpr double RTOL = 1e-2; - constexpr double ATOL = 1e-8; + constexpr double RTOL = 1e-5; + constexpr double ATOL = 1e-5; bool distClose(double a, double b) { @@ -61,6 +62,7 @@ class EquidepthTest : public ::testing::Test std::vector data, query; std::vector> gtDist; + std::vector gtIds; void SetUp() override { @@ -73,6 +75,7 @@ class EquidepthTest : public ::testing::Test std::vector gtI((size_t)NQ * K); std::vector gtD((size_t)NQ * K); bf.searchIndex(query.data(), NQ, K, gtI.data(), gtD.data()); + gtIds = gtI; gtDist.assign(NQ, {}); for (int q = 0; q < NQ; q++) { @@ -100,18 +103,48 @@ class EquidepthTest : public ::testing::Test ours[j] = {D[(size_t)q * K + j], I[(size_t)q * K + j]}; std::sort(ours.begin(), ours.end()); - int uniq = 0; + std::set returned; for (int j = 0; j < K; j++) - if (j == 0 || ours[j].second != ours[j - 1].second || - !distClose(ours[j].first, ours[j - 1].first)) - uniq++; - - int compareCount = std::min(uniq, K); - for (int j = 0; j < compareCount; j++) + { + ASSERT_LT(ours[j].second, static_cast(N)); + EXPECT_TRUE(returned.insert(ours[j].second).second) + << "duplicate ID, query " << q << ", rank " << j; + double actual = 0.0; + for (int d = 0; d < DIM; ++d) + { + double delta = static_cast(data[ours[j].second * DIM + d]) - query[q * DIM + d]; + actual += delta * delta; + } + EXPECT_TRUE(distClose(ours[j].first, actual)) + << "returned ID/distance mismatch, query " << q << ", rank " << j; EXPECT_TRUE(distClose(ours[j].first, gtDist[q][j])) << "query " << q << " neighbour " << j << ": got " << ours[j].first << ", exact GT " << gtDist[q][j] << " (equi-depth kNN not exact on non-normalized data)"; + } + for (int j = 0; j < K; ++j) + EXPECT_EQ(returned.count(gtIds[q * K + j]), 1u) + << "missing exact neighbor, query " << q << ", rank " << j; + } + } + + void checkSearchRunsOnNonNormalized(daisy::SimilaritySearchAlgorithm *algo) + { + algo->setNormalized(false); + EXPECT_FALSE(algo->getNormalized()) << "setNormalized(false) did not take effect"; + + algo->buildIndex(data.data(), N, DIM); + algo->setNumThreads(1); + + std::vector I((size_t)NQ * K); + std::vector D((size_t)NQ * K); + algo->searchIndex(query.data(), NQ, K, I.data(), D.data()); + + for (size_t i = 0; i < D.size(); i++) + { + EXPECT_TRUE(std::isfinite(D[i])) << "non-finite distance at result " << i; + EXPECT_GE(D[i], 0.0f) << "negative distance at result " << i; + EXPECT_LT(I[i], static_cast(N)) << "out-of-range index at result " << i; } } }; @@ -131,7 +164,7 @@ TEST_F(EquidepthTest, Messi) TEST_F(EquidepthTest, Fresh) { daisy::Fresh algo(daisy::DistanceType::L2_SQUARED); - checkExactOnNonNormalized(&algo); + checkSearchRunsOnNonNormalized(&algo); } // Hercules only supports z-normalized data: setNormalized(false) must be rejected. @@ -147,10 +180,8 @@ TEST_F(EquidepthTest, HerculesRejectsNonNormalized) // index's lower bounds were only correct while no other index had installed its own. The // query paths of Messi and LbBruteforce now pass their own tables down explicitly. // -// These tests assert exactly what that guarantees: an index answers a query identically -// whether or not a second index with different breakpoints exists. They deliberately do not -// compare against brute force -- MESSI's equi-depth kNN is not exact on every dataset, and -// that is a separate matter from whether a second index perturbs it. +// Each index must remain exact against brute force whether or not a second +// index with different breakpoints exists. class TwoIndexEquidepthTest : public ::testing::Test { protected: @@ -174,6 +205,12 @@ class TwoIndexEquidepthTest : public ::testing::Test template void expectSecondIndexDoesNotPerturbFirst() { + daisy::BruteForceSearch brute(daisy::DistanceType::L2_SQUARED); + brute.buildIndex(dataA.data(), N, DIM); + std::vector expected_I((size_t)NQ * K); + std::vector expected_D((size_t)NQ * K); + brute.searchIndex(queryA.data(), NQ, K, expected_I.data(), expected_D.data()); + Algo a(daisy::DistanceType::L2_SQUARED); a.setNormalized(false); a.buildIndex(dataA.data(), N, DIM); @@ -200,13 +237,28 @@ class TwoIndexEquidepthTest : public ::testing::Test EXPECT_FLOAT_EQ(after_D[i], alone_D[i]) << "result " << i << " changed once a second index existed"; } + for (int q = 0; q < NQ; ++q) + { + std::set alone(alone_I.begin() + q * K, alone_I.begin() + (q + 1) * K); + std::set after(after_I.begin() + q * K, after_I.begin() + (q + 1) * K); + EXPECT_EQ(alone.size(), static_cast(K)); + EXPECT_EQ(after.size(), static_cast(K)); + for (int j = 0; j < K; ++j) + { + EXPECT_EQ(alone.count(expected_I[q * K + j]), 1u); + EXPECT_EQ(after.count(expected_I[q * K + j]), 1u); + EXPECT_TRUE(distClose(alone_D[q * K + j], expected_D[q * K + j])); + EXPECT_TRUE(distClose(after_D[q * K + j], expected_D[q * K + j])); + } + } } }; -// Messi is deliberately not covered here: in equi-depth mode its kNN output already -// drifts between two identical consecutive queries on a single index (measured at -// ~9/200 results, unchanged by this migration and independent of worker count), so an -// exact-equality assertion would be flaky for reasons unrelated to breakpoint plumbing. +TEST_F(TwoIndexEquidepthTest, MessiQueriesAreUnaffectedByASecondIndex) +{ + expectSecondIndexDoesNotPerturbFirst(); +} + TEST_F(TwoIndexEquidepthTest, LbBruteforceQueriesAreUnaffectedByASecondIndex) { expectSecondIndexDoesNotPerturbFirst(); diff --git a/tests/test_messi_query_layout.py b/tests/test_messi_query_layout.py new file mode 100644 index 0000000..f230275 --- /dev/null +++ b/tests/test_messi_query_layout.py @@ -0,0 +1,88 @@ +"""Exact MESSI queries across NumPy layouts, including header-strided fvecs views.""" +import unittest +import numpy as np +import daisy + + +class MessiQueryLayoutTest(unittest.TestCase): + def test_layouts_static_and_streaming(self): + rng = np.random.default_rng(73) + # Sparse positive descriptors expose the lowest adaptive SAX interval. + database = rng.integers(0, 80, size=(1024, 128)).astype(np.float32) + database[rng.random(database.shape) < 0.7] = 0 + queries = rng.integers(0, 80, size=(8, 128)).astype(np.float32) + queries[rng.random(queries.shape) < 0.7] = 0 + header_rows = np.empty((8, 129), dtype=np.float32) + header_rows[:, 0] = 128 + header_rows[:, 1:] = queries + alternate_rows = np.empty((16, 128), dtype=np.float32) + alternate_rows[::2] = queries + layouts = { + "dense": queries, + "fvecs": header_rows[:, 1:], + "row_slice": alternate_rows[::2], + "fortran": np.asfortranarray(queries), + "reversed_rows": queries[::-1], + "reversed_columns": queries[:, ::-1], + } + for normalized in (False, True): + for workers in (1, 4): + index = daisy.Messi(daisy.DistanceType.L2_SQUARED) + index.setNormalized(normalized) + index.setIndexWorkers(2) + index.setSearchWorkers(workers) + initial = np.ascontiguousarray(database[:768]) + index.buildIndex(initial) + for size in (768, 1024): + if size == 1024: + index.insertBatch(database[768:]) + for name, query in layouts.items(): + exact = ((database[:size].astype(np.float64)[None, :, :] - + query.astype(np.float64)[:, None, :]) ** 2).sum(axis=2) + expected = np.sort(exact, axis=1)[:, :10] + config = daisy.SearchConfig() + config.type = daisy.QueryType.TOP_K + config.k = 10 + for api in (10, config): + with self.subTest(normalized=normalized, workers=workers, + size=size, layout=name, api=str(api)): + ids, distances = index.searchIndex(query, api) + ids = np.asarray(ids, dtype=np.int64) + distances = np.asarray(distances) + np.testing.assert_array_equal(distances, expected) + np.testing.assert_array_equal( + distances, np.take_along_axis(exact, ids, axis=1)) + for row in ids: + self.assertEqual(len(set(row)), 10) + radius = float(np.median(expected[:, -1])) + config.type = daisy.QueryType.RANGE + config.r = radius + ids, distances = index.searchIndex(query, config) + for q, (row_ids, row_d) in enumerate(zip(ids, distances)): + self.assertEqual(set(row_ids), set(np.flatnonzero(exact[q] <= radius))) + np.testing.assert_array_equal(row_d, exact[q, row_ids]) + self.assertEqual(len(row_ids), len(set(row_ids))) + + def test_invalid_query_shape_and_empty_index(self): + index = daisy.Messi(daisy.DistanceType.L2_SQUARED) + config = daisy.SearchConfig() + config.k = 1 + for api in (1, config): + with self.assertRaisesRegex(RuntimeError, "built"): + index.searchIndex(np.zeros((1, 128), np.float32), api) + database = np.arange(32 * 128, dtype=np.float32).reshape(32, 128) + index.setNormalized(False) + index.setIndexWorkers(1) + index.buildIndex(database) + for api in (1, config): + with self.assertRaisesRegex(RuntimeError, "dimension"): + index.searchIndex(np.zeros((2, 127), np.float32), api) + with self.assertRaisesRegex(RuntimeError, "2D"): + index.searchIndex(np.zeros(128, np.float32), api) + config.k = 0 + with self.assertRaisesRegex(RuntimeError, "positive"): + index.searchIndex(np.zeros((1, 128), np.float32), config) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_sax_breakpoints.cpp b/tests/test_sax_breakpoints.cpp new file mode 100644 index 0000000..b9c0572 --- /dev/null +++ b/tests/test_sax_breakpoints.cpp @@ -0,0 +1,68 @@ +#include +#include +#include +#include +#include + +namespace { +// Restore the process-wide table even if a fatal assertion exits the test. +struct ActiveTableGuard { + const float *cuts = daisy::daisy_active_breakpoints; + const float *max_cuts = daisy::daisy_active_breakpoints_max; + ~ActiveTableGuard() { daisy::set_active_breakpoints(cuts, max_cuts); } +}; +} + +TEST(SaxBreakpointTest, ClassifiesShiftedRepeatedAndBoundaryCuts) +{ + ActiveTableGuard guard; + for (int cardinality : {2, 4, 256}) { + const int offset = (cardinality - 1) * (cardinality - 2) / 2; + for (float shift : {-100.0f, 0.0f, 100.0f}) { + // Poison the previous cardinality's last cut: it is not a lower + // endpoint of this cardinality's first interval. + std::vector table(offset + cardinality - 1, 10000.0f); + for (int i = 0; i < cardinality - 1; ++i) + table[offset + i] = shift + i / 2; + daisy::set_active_breakpoints(table.data(), nullptr); + std::vector probes{shift - 1, shift, shift + 0.25f, + table.back(), table.back() + 1}; + for (int i = 0; i < cardinality - 1; ++i) + probes.push_back(table[offset + i]); + for (float value : probes) { + int expected = 0; + while (expected < cardinality - 1 && value > table[offset + expected]) + ++expected; + daisy::sax_type from_paa = 0, from_ts = 0; + std::array series; + series.fill(value); + ASSERT_EQ(daisy::sax_from_paa(&value, &from_paa, 1, cardinality, 8), daisy::SUCCESS); + ASSERT_EQ(daisy::sax_from_ts(series.data(), &from_ts, 8, 1, cardinality, 8), daisy::SUCCESS); + EXPECT_EQ(from_paa, expected) << "cardinality=" << cardinality << " value=" << value; + EXPECT_EQ(from_ts, expected) << "cardinality=" << cardinality << " value=" << value; + } + } + } +} + +TEST(SaxBreakpointTest, LowPositivePaaNeverOverestimatesRecordDistance) +{ + ActiveTableGuard guard; + constexpr int cardinality = 256, segments = 16, values_per_segment = 8; + constexpr int offset = (cardinality - 1) * (cardinality - 2) / 2; + std::vector table(offset + cardinality - 1, 200.0f); + for (int i = 0; i < cardinality - 1; ++i) + table[offset + i] = 2.0f + i / 4; + daisy::set_active_breakpoints(table.data(), nullptr); + std::array record_paa, query_paa; + record_paa.fill(1.5f); + query_paa.fill(3.0f); + std::array sax, bits; + bits.fill(8); + ASSERT_EQ(daisy::sax_from_paa(record_paa.data(), sax.data(), segments, cardinality, 8), daisy::SUCCESS); + const float exact = segments * values_per_segment * (3.0f - 1.5f) * (3.0f - 1.5f); + EXPECT_LE(daisy::minidist_paa_to_isax(query_paa.data(), sax.data(), bits.data(), + 8, cardinality, segments, -2000000, 2000000, values_per_segment, table.data()), exact); + EXPECT_LE(daisy::minidist_paa_to_isax_rawa_SIMD(query_paa.data(), sax.data(), bits.data(), + 8, cardinality, segments, -2000000, 2000000, values_per_segment, table.data()), exact); +}