From 0fd3538ad98fec2a31990de060dc5cc96e3eee27 Mon Sep 17 00:00:00 2001 From: sophisid Date: Mon, 14 Sep 2026 22:06:33 +0200 Subject: [PATCH 1/4] streaming (incremental) MESSI --- CONTINUITY.md | 5 +- README.md | 19 +- benchmark/CMakeLists.txt | 12 ++ benchmark/bm_Messi_Streaming.cpp | 59 +++++++ demos/CMakeLists.txt | 7 + demos/demo_Messi_Streaming.cpp | 42 +++++ demos/demo_Messi_Streaming.py | 23 +++ docs/demos-guide.md | 5 +- docs/how-to-contribute.md | 2 +- lib/algos/Messi.cpp | 224 +++++++++++++++++++++++- lib/algos/Messi.hpp | 23 ++- lib/algos/SimilaritySearchAlgorithm.hpp | 2 +- lib/isax/iSAXIndex.cpp | 114 +++++++----- lib/isax/iSAXIndex.hpp | 11 +- pybinds/setup.cpp | 36 +++- tests/CMakeLists.txt | 24 +++ tests/test_Messi_Streaming.cpp | 214 ++++++++++++++++++++++ tests/test_streaming.py | 60 ++++++- 18 files changed, 812 insertions(+), 70 deletions(-) create mode 100644 benchmark/bm_Messi_Streaming.cpp create mode 100644 demos/demo_Messi_Streaming.cpp create mode 100644 demos/demo_Messi_Streaming.py create mode 100644 tests/test_Messi_Streaming.cpp diff --git a/CONTINUITY.md b/CONTINUITY.md index c67f0b7..0ba0bd5 100644 --- a/CONTINUITY.md +++ b/CONTINUITY.md @@ -65,10 +65,13 @@ These algorithms do not all share the same dependency profile or feature set. - Top-k search is the baseline capability. - Range search is modeled through `SearchConfig` and is not universally implemented. -- Streaming insert is implemented by `BruteForceSearch`, `LbBruteforce`, and `Coconut`; +- Streaming insert is implemented by `BruteForceSearch`, `LbBruteforce`, `Messi`, and `Coconut`; the base implementation still throws for algorithms that do not support it. - Bruteforce streaming grows the owned in-memory database incrementally. LbBruteforce also computes SAX summaries incrementally, using the breakpoint set fixed by the initial build. +- MESSI uses those fixed breakpoints to route stable SAX/position records into its live iSAX + tree, including new-root creation and leaf splitting. Callers must serialize inserts and + queries; simultaneous update/search is not supported. - `setNormalized(bool)` is a declaration about the input data, not a preprocessing step. - Data can come from in-memory arrays or file-backed sources via `DataSource`. diff --git a/README.md b/README.md index 124d7cb..b7e012d 100644 --- a/README.md +++ b/README.md @@ -48,7 +48,7 @@ The following table summarizes the key features of each algorithm: |-----------|-------------| | **Bruteforce** | Naive parallel similarity search implementation with incremental streaming inserts | | **Lower Bound Bruteforce** | Optimized bruteforce with lower bounding and incremental streaming inserts | -| **[MESSI](https://helios2.mi.parisdescartes.fr/~themisp/messi/)** | In-memory parallel similarity search | +| **[MESSI](https://helios2.mi.parisdescartes.fr/~themisp/messi/)** | In-memory parallel similarity search with incremental iSAX-tree inserts | | **[PARIS](https://helios2.mi.parisdescartes.fr/~themisp/paris/)** | Disk-based parallel similarity search | | **[SING](https://helios2.mi.parisdescartes.fr/~themisp/sing/)** | GPU-accelerated in-memory parallel similarity search | | **[Odyssey](https://helios2.mi.parisdescartes.fr/~themisp/odyssey/)** | Distributed and parallel in-memory similarity search | @@ -60,20 +60,24 @@ The following table summarizes the key features of each algorithm: ### Incremental streaming inserts -`BruteForceSearch`, `LbBruteforce`, and `Coconut` implement the common streaming API. Build +`BruteForceSearch`, `LbBruteforce`, `Messi`, and `Coconut` implement the common streaming API. Build the initial index once, then append one series or a contiguous batch without rebuilding: ```cpp -daisy::BruteForceSearch search(daisy::DistanceType::L2_SQUARED); +daisy::Messi search(daisy::DistanceType::L2_SQUARED); search.buildIndex(initial_data, initial_size, dim); search.insert(one_series); search.insertBatch(batch_data, batch_size); ``` Inserted series receive consecutive IDs beginning at the size of the initial database and are -immediately visible to top-k and range searches. `LbBruteforce` computes a SAX summary for each -insert using the breakpoints established during the initial build. Inserts can reallocate the -owned database, so callers should not retain a pointer returned by `getDatabase()` across them. +immediately visible to supported top-k and range searches. `LbBruteforce` and `Messi` compute a +SAX summary for each insert using the breakpoints established during the initial build. MESSI +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. + +Streaming updates are not concurrent with queries. Inserts can reallocate the owned database, +so callers should not retain a pointer returned by `getDatabase()` across them. @@ -195,6 +199,7 @@ cd build ./benchmark/bm_bruteforce_L2Square ./benchmark/bm_LbBruteforce_L2Square ./benchmark/bm_Messi_L2Square +./benchmark/bm_Messi_Streaming # Advanced algorithms (if available) ./benchmark/bm_Odyssey_L2Square # MPI required @@ -221,5 +226,3 @@ For questions and suggestions through mail, you can contact us at [manos.chatzak - - diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt index 05f951f..aa8a5f0 100644 --- a/benchmark/CMakeLists.txt +++ b/benchmark/CMakeLists.txt @@ -209,6 +209,18 @@ if(DEBUG_MSG) message(STATUS "Include directories added for bm_Messi_L2Square.") endif() +add_executable(bm_Messi_Streaming bm_Messi_Streaming.cpp) +target_link_libraries(bm_Messi_Streaming + PRIVATE + benchmark::benchmark + benchmark::benchmark_main + dino_lib +) +target_include_directories(bm_Messi_Streaming + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib +) + # ////// FAISS FLAT ////// if(DEBUG_MSG) message(STATUS "---") diff --git a/benchmark/bm_Messi_Streaming.cpp b/benchmark/bm_Messi_Streaming.cpp new file mode 100644 index 0000000..8aaebdc --- /dev/null +++ b/benchmark/bm_Messi_Streaming.cpp @@ -0,0 +1,59 @@ +#include + +#include "../lib/algos/Messi.hpp" + +#include +#include + +namespace +{ + constexpr int DIM = 96; + constexpr int INITIAL = 4096; + + std::vector makeSeries(int n, int phase_offset) + { + std::vector data(static_cast(n) * DIM); + for (int i = 0; i < n; ++i) + { + float *series = data.data() + static_cast(i) * DIM; + const float phase = static_cast(phase_offset + i) * 0.013f; + for (int j = 0; j < DIM; ++j) + series[j] = std::sin(0.07f * j + phase) + + 0.5f * std::cos(0.19f * j - phase); + } + return data; + } +} + +static void BM_Messi_InsertBatch(benchmark::State &state) +{ + state.PauseTiming(); + const daisy::idx_t batch_size = static_cast(state.range(0)); + auto initial = makeSeries(INITIAL, 0); + auto batch = makeSeries(static_cast(batch_size), INITIAL); + + daisy::MessiConfig config; + config.index_workers = 2; + config.search_workers = 2; + daisy::Messi search(daisy::DistanceType::L2_SQUARED, config); + search.buildIndex(initial.data(), INITIAL, DIM); + state.ResumeTiming(); + + for (auto _ : state) + { + search.insertBatch(batch.data(), batch_size); + benchmark::ClobberMemory(); + } + + state.SetItemsProcessed( + static_cast(state.iterations()) * static_cast(batch_size)); +} + +BENCHMARK(BM_Messi_InsertBatch) + ->Arg(1) + ->Arg(64) + ->Arg(1024) + ->Iterations(20) + ->Unit(benchmark::kMicrosecond); + +BENCHMARK_MAIN(); diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt index 6778c9d..209426c 100644 --- a/demos/CMakeLists.txt +++ b/demos/CMakeLists.txt @@ -204,6 +204,13 @@ if(BUILD_DEMO) message(STATUS "Include directories added for demo_Messi_L2Square.") endif() + add_executable(demo_Messi_Streaming demo_Messi_Streaming.cpp) + target_link_libraries(demo_Messi_Streaming PRIVATE dino_lib commons_lib) + target_include_directories(demo_Messi_Streaming PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + # ////// MESSI DTW ////// if(DEBUG_MSG) message(STATUS "---") diff --git a/demos/demo_Messi_Streaming.cpp b/demos/demo_Messi_Streaming.cpp new file mode 100644 index 0000000..40fc336 --- /dev/null +++ b/demos/demo_Messi_Streaming.cpp @@ -0,0 +1,42 @@ +// MESSI streaming: build once, then update the live iSAX tree without rebuilding. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include + +int main() +{ + const daisy::idx_t dim = 96; + const daisy::idx_t initial = 5000; + const daisy::idx_t batch = 1000; + const daisy::idx_t n_query = 5; + const daisy::idx_t k = 5; + + float *stream = loadRandomData(initial + batch, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + daisy::Messi search(daisy::DistanceType::L2_SQUARED); + search.setIndexWorkers(2); + search.setSearchWorkers(4); + search.buildIndex(stream, initial, dim); + + search.insert(stream + initial * dim); + search.insertBatch(stream + (initial + 1) * dim, batch - 1); + + daisy::idx_t *indices = new daisy::idx_t[n_query * k]; + float *distances = new float[n_query * k]; + search.searchIndex(query, n_query, k, indices, distances); + + std::printf("MESSI now contains %llu series. Query 0 kNN: ", + search.getNDatabase()); + for (daisy::idx_t j = 0; j < k; ++j) + std::printf("%llu(%.3f) ", indices[j], distances[j]); + std::printf("\n"); + + delete[] stream; + delete[] query; + delete[] indices; + delete[] distances; + return 0; +} diff --git a/demos/demo_Messi_Streaming.py b/demos/demo_Messi_Streaming.py new file mode 100644 index 0000000..e416160 --- /dev/null +++ b/demos/demo_Messi_Streaming.py @@ -0,0 +1,23 @@ +"""Build MESSI once and append streaming batches without rebuilding.""" + +import numpy as np + +from daisy import DistanceType, Messi + + +rng = np.random.default_rng(7) +stream = rng.normal(size=(6000, 96)).astype(np.float32) +queries = rng.normal(size=(5, 96)).astype(np.float32) + +index = Messi(DistanceType.L2_SQUARED) +index.setIndexWorkers(2) +index.setSearchWorkers(4) +index.buildIndex(stream[:5000]) + +index.insert(stream[5000]) +index.insertBatch(stream[5001:]) + +indices, distances = index.searchIndex(queries, 5) +print("MESSI size:", 6000) +print("Query 0 IDs:", indices[0]) +print("Query 0 distances:", distances[0]) diff --git a/docs/demos-guide.md b/docs/demos-guide.md index f825f2a..6b80be5 100644 --- a/docs/demos-guide.md +++ b/docs/demos-guide.md @@ -4,9 +4,10 @@ The demos module provides practical examples of how to use the DaiSy library's a Each demo illustrates a specific algorithm or specific distance metric. This module includes both C++ and Python implementations for various algorithms and use cases. Most demos follow the same batch pattern: `buildIndex(...)` once, then `searchIndex(...)`. -**Bruteforce**, **LbBruteforce**, and **Coconut** additionally support streaming through +**Bruteforce**, **LbBruteforce**, **MESSI**, and **Coconut** additionally support streaming through `insert(...)` and `insertBatch(...)`. See `demo_Bruteforce_Streaming`, -`demo_LbBruteforce_Streaming`, and `demo_Coconut_Streaming` for live-index examples. +`demo_LbBruteforce_Streaming`, `demo_Messi_Streaming`, and `demo_Coconut_Streaming` for +live-index examples. ## Demo Program Structure diff --git a/docs/how-to-contribute.md b/docs/how-to-contribute.md index adc1939..6794b9e 100644 --- a/docs/how-to-contribute.md +++ b/docs/how-to-contribute.md @@ -12,7 +12,7 @@ Here is a (non-exhaustive) mockup of our future and ongoing goals: - Extension of DaiSy for subsequence similarity search - Extension of DaiSy for more algorithms (e.g., SFA, Hercules, Dumpy, etc.) -- Streaming / updatable indexing for more algorithms (currently supported by Bruteforce, LbBruteforce, and Coconut) +- Streaming / updatable indexing for more algorithms (currently supported by Bruteforce, LbBruteforce, MESSI, and Coconut) - Implementation of a DaiSy autotuner to automatically optimize indexing and search parameters - Extension of DaiSy to support learned optimization approaches, e.g., LeaFi and ProS diff --git a/lib/algos/Messi.cpp b/lib/algos/Messi.cpp index 2eb6fcf..345ba68 100644 --- a/lib/algos/Messi.cpp +++ b/lib/algos/Messi.cpp @@ -3,6 +3,11 @@ #include #include #include +#include +#include +#include +#include +#include #include #include @@ -637,6 +642,75 @@ namespace daisy this->paa_segments = config.paa_segments; } + void Messi::reserveDatabase(idx_t required_capacity) + { + if (required_capacity <= this->database_capacity && this->owns_database) + return; + + idx_t new_capacity = std::max(this->database_capacity, 1); + while (new_capacity < required_capacity) + { + if (new_capacity > std::numeric_limits::max() / 2) + { + new_capacity = required_capacity; + break; + } + new_capacity *= 2; + } + + if (this->dim == 0 || + new_capacity > std::numeric_limits::max() / this->dim) + throw std::length_error("Messi database is too large"); + + std::unique_ptr grown_database( + new float[static_cast(new_capacity) * + static_cast(this->dim)]); + if (this->n_database > 0) + { + std::copy_n(this->database, + static_cast(this->n_database) * + static_cast(this->dim), + grown_database.get()); + } + + if (this->owns_database) + delete[] this->database; + this->database = grown_database.release(); + this->owns_database = true; + this->database_capacity = new_capacity; + } + + void Messi::reserveSaxCache(idx_t required_capacity) + { + if (required_capacity <= this->sax_cache_capacity) + return; + + idx_t new_capacity = std::max(this->sax_cache_capacity, 1); + while (new_capacity < required_capacity) + { + if (new_capacity > std::numeric_limits::max() / 2) + { + new_capacity = required_capacity; + break; + } + new_capacity *= 2; + } + + const size_t segments = static_cast(this->index->settings->paa_segments); + if (segments == 0 || + new_capacity > std::numeric_limits::max() / segments) + throw std::length_error("Messi SAX cache is too large"); + + void *grown = std::realloc( + this->index->sax_cache, + static_cast(new_capacity) * segments * sizeof(sax_type)); + if (grown == nullptr) + throw std::bad_alloc(); + + this->index->sax_cache = static_cast(grown); + this->sax_cache_capacity = new_capacity; + } + void *indexCreationWorker(void *transferdata) { sax_type *sax = (sax_type *)malloc(sizeof(sax_type) * ((buffer_data_inmemory *)transferdata)->index->settings->paa_segments); @@ -740,9 +814,19 @@ namespace daisy void Messi::buildIndex(DataSource *data_source) { + if (data_source == nullptr) + throw std::invalid_argument("Messi::buildIndex received a null data source"); + if (this->index != nullptr) + throw std::runtime_error("Messi::buildIndex may only be called once per instance"); + this->dim = data_source->getDim(); this->n_database = data_source->getTotalRecords(); + if (this->dim == 0) + throw std::invalid_argument("Messi::buildIndex requires a positive dimension"); + if (this->index_workers < 1) + throw std::invalid_argument("Messi::buildIndex requires at least one index worker"); + if (this->n_database == 0) { @@ -780,6 +864,7 @@ namespace daisy delete[] record; this->owns_database = true; } + this->database_capacity = this->n_database; this->index_settings = isax_index_settings_init("", this->dim, @@ -811,7 +896,12 @@ namespace daisy int node_counter = 0; pthread_t threadid[this->index_workers]; buffer_data_inmemory *input_data = (buffer_data_inmemory *)malloc(sizeof(buffer_data_inmemory) * (this->index_workers)); - index->sax_cache = (sax_type *)malloc(sizeof(sax_type) * index->settings->paa_segments * this->n_database); + this->sax_cache_capacity = std::max(this->n_database, 1); + index->sax_cache = (sax_type *)malloc( + sizeof(sax_type) * index->settings->paa_segments * + static_cast(this->sax_cache_capacity)); + if (index->sax_cache == nullptr) + throw std::bad_alloc(); pthread_barrier_t lock_barrier1, lock_barrier2; pthread_barrier_init(&lock_barrier1, NULL, this->index_workers + 1); pthread_barrier_init(&lock_barrier2, NULL, this->index_workers + 1); @@ -870,6 +960,132 @@ namespace daisy free(nodesize); } + void Messi::insert(const float *series) + { + insertBatch(series, 1); + } + + void Messi::insertBatch(const float *data, idx_t n) + { + if (n == 0) + return; + if (this->index == nullptr || this->index_settings == nullptr || this->dim == 0) + throw std::runtime_error("Messi::insertBatch requires an initial buildIndex first"); + if (data == nullptr) + throw std::invalid_argument("Messi::insertBatch received null data"); + if (n > std::numeric_limits::max() - this->n_database) + throw std::length_error("Messi database size overflow"); + + const idx_t required_capacity = this->n_database + n; + if (required_capacity > std::numeric_limits::max() / this->dim) + throw std::length_error("Messi position offset overflow"); + if (required_capacity > std::numeric_limits::max()) + throw std::length_error("Messi SAX cache size overflow"); + if (n > std::numeric_limits::max() / this->dim) + throw std::length_error("Messi insert batch is too large"); + + const size_t current_values = + static_cast(this->n_database) * static_cast(this->dim); + const size_t inserted_values = + static_cast(n) * static_cast(this->dim); + const uintptr_t database_begin = reinterpret_cast(this->database); + const uintptr_t database_end = + database_begin + current_values * sizeof(float); + const uintptr_t data_address = reinterpret_cast(data); + const bool aliases_database = + current_values > 0 && data_address >= database_begin && data_address < database_end; + size_t source_offset = 0; + if (aliases_database) + { + const uintptr_t byte_offset = data_address - database_begin; + if (byte_offset % sizeof(float) != 0) + throw std::invalid_argument("Messi::insertBatch received an unaligned database pointer"); + source_offset = static_cast(byte_offset / sizeof(float)); + if (inserted_values > current_values - source_offset) + throw std::invalid_argument("Messi::insertBatch source exceeds the live database"); + } + + // MESSI may initially borrow a caller-owned contiguous buffer. The first + // update converts that view to owned, growable storage without rebuilding + // the tree; existing position offsets remain valid after the copy. + reserveDatabase(required_capacity); + reserveSaxCache(required_capacity); + const float *source_data = aliases_database ? this->database + source_offset : data; + + const size_t segments = static_cast(this->index->settings->paa_segments); + if (n > std::numeric_limits::max() / segments) + throw std::length_error("Messi SAX insert batch is too large"); + + IncrementalRecordBlock block; + block.sax = std::make_unique(static_cast(n) * segments); + block.positions = std::make_unique(static_cast(n)); + + activateBreakpoints(); + for (idx_t i = 0; i < n; ++i) + { + sax_type *sax = block.sax.get() + static_cast(i) * segments; + const float *series = source_data + static_cast(i) * this->dim; + if (!this->distance_computer->compute_sax_from_ts( + series, + sax, + this->index->settings->ts_values_per_paa_segment, + this->index->settings->paa_segments, + this->index->settings->sax_alphabet_cardinality, + this->index->settings->sax_bit_cardinality)) + throw std::runtime_error("Messi::insertBatch failed to compute SAX representation"); + + block.positions[static_cast(i)] = + static_cast(this->n_database + i) * this->dim; + } + + // Resolve every destination root before inserting any record. New roots + // are registered in the FBL lookup used by MESSI's approximate search, + // but the records themselves stay in the stable block above. + std::vector roots(static_cast(n)); + pthread_mutex_t lock_firstnode = PTHREAD_MUTEX_INITIALIZER; + auto *fbl = reinterpret_cast(this->index->fbl); + for (idx_t i = 0; i < n; ++i) + { + sax_type *sax = block.sax.get() + static_cast(i) * segments; + root_mask_type root_mask = 0; + CREATE_MASK(root_mask, this->index, sax); + roots[static_cast(i)] = get_or_create_pRecBuf_root( + fbl, root_mask, this->index, &lock_firstnode, this->index_workers); + if (roots[static_cast(i)] == nullptr) + { + pthread_mutex_destroy(&lock_firstnode); + throw std::bad_alloc(); + } + } + pthread_mutex_destroy(&lock_firstnode); + + this->incremental_record_blocks.reserve( + this->incremental_record_blocks.size() + 1); + this->incremental_record_blocks.push_back(std::move(block)); + IncrementalRecordBlock &stored = this->incremental_record_blocks.back(); + + std::copy_n(source_data, + inserted_values, + this->database + static_cast(this->n_database) * this->dim); + std::copy_n(stored.sax.get(), + static_cast(n) * segments, + this->index->sax_cache + static_cast(this->n_database) * segments); + + for (idx_t i = 0; i < n; ++i) + { + isax_node_record record{}; + record.sax = stored.sax.get() + static_cast(i) * segments; + record.position = stored.positions.get() + static_cast(i); + record.ts = nullptr; + record.insertion_mode = static_cast(NO_TMP | PARTIAL); + add_record_to_node(this->index, roots[static_cast(i)], &record, 1); + } + + this->n_database = required_capacity; + this->index->total_records += n; + this->index->sax_cache_size = static_cast(required_capacity); + } + void Messi::searchIndexL2Squared(const float *query, const idx_t n_query, const idx_t k, idx_t *I, float *D) { ts_type *paa = (ts_type *)malloc(sizeof(ts_type) * index->settings->paa_segments); @@ -1079,6 +1295,7 @@ namespace daisy std::vector> &I, std::vector> &D) { + activateBreakpoints(); if (config.type == QueryType::TOP_K) { SimilaritySearchAlgorithm::searchIndex(query, n_query, config, I, D); @@ -1122,6 +1339,7 @@ namespace daisy void Messi::searchIndex(const float *query, const idx_t n_query, const idx_t k, idx_t *I, float *D) { + activateBreakpoints(); if (this->distance_type == DistanceType::L2_SQUARED) { searchIndexL2Squared(query, n_query, k, I, D); @@ -1358,6 +1576,8 @@ namespace daisy if (index_settings != nullptr) { + if (daisy_active_breakpoints == index_settings->breakpoints) + set_active_breakpoints(nullptr, nullptr); if (index_settings->bit_masks != nullptr) { free(index_settings->bit_masks); @@ -1366,6 +1586,8 @@ namespace daisy { free(index_settings->max_sax_cardinalities); } + free(index_settings->breakpoints_owned); + free(index_settings->breakpoints_max_owned); free(index_settings); } } diff --git a/lib/algos/Messi.hpp b/lib/algos/Messi.hpp index 9cf8034..4fac58b 100644 --- a/lib/algos/Messi.hpp +++ b/lib/algos/Messi.hpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -79,6 +80,21 @@ namespace daisy int index_workers = 2; int n_pqueue = 42; bool owns_database = false; // track ownership of database buffer + idx_t database_capacity = 0; + idx_t sax_cache_capacity = 0; + + // Tree leaf buffers keep raw pointers to SAX words and positions. Each + // incremental batch therefore owns stable backing arrays for as long as + // the MESSI instance is alive; growing the vector only moves unique_ptrs. + struct IncrementalRecordBlock + { + std::unique_ptr sax; + std::unique_ptr positions; + }; + std::vector incremental_record_blocks; + + void reserveDatabase(idx_t required_capacity); + void reserveSaxCache(idx_t required_capacity); pqueue_bsf MESSI_search_topk_L2Squared(ts_type *ts, ts_type *paa, node_list *nodelist, idx_t k); pqueue_bsf MESSI_search_topk_DTW(ts_type *ts, node_list *nodelist, idx_t k); @@ -103,6 +119,11 @@ namespace daisy throw std::runtime_error("Messi requires in-memory data. Use buildIndex(database, n_database, dim) instead."); } + // Callers must not overlap incremental updates with searches. The + // breakpoints established by buildIndex remain fixed for all inserts. + void insert(const float *series) override; + void insertBatch(const float *data, idx_t n) override; + void searchIndex(const float *query, const idx_t n_query, const idx_t k, idx_t *I, float *D) override; void searchIndex(const float *query, idx_t n_query, const SearchConfig &config, @@ -140,7 +161,7 @@ namespace daisy void setReadBlockLength(int n) { read_block_length = n; } int getWarpingWindow() const { return warping_window; } - ~Messi(); + ~Messi() override; }; } diff --git a/lib/algos/SimilaritySearchAlgorithm.hpp b/lib/algos/SimilaritySearchAlgorithm.hpp index defe691..e982708 100644 --- a/lib/algos/SimilaritySearchAlgorithm.hpp +++ b/lib/algos/SimilaritySearchAlgorithm.hpp @@ -104,7 +104,7 @@ namespace daisy } // Streaming API: incrementally add series to a live index. Only some algorithms - // support it (BruteForceSearch, LbBruteforce, and Coconut); the default throws. + // support it (BruteForceSearch, LbBruteforce, Messi, and Coconut); the default throws. virtual void insert(const float *series) { (void)series; diff --git a/lib/isax/iSAXIndex.cpp b/lib/isax/iSAXIndex.cpp index a23a152..02dd685 100644 --- a/lib/isax/iSAXIndex.cpp +++ b/lib/isax/iSAXIndex.cpp @@ -355,6 +355,10 @@ namespace daisy { fbl->soft_buffers[i].initialized = 0; fbl->soft_buffers[i].finished = 0; + fbl->soft_buffers[i].node = NULL; + fbl->soft_buffers[i].sax_records = NULL; + fbl->soft_buffers[i].pos_records = NULL; + fbl->soft_buffers[i].max_buffer_size = NULL; fbl->soft_buffers[i].buffer_size = NULL; } return fbl; @@ -572,6 +576,68 @@ namespace daisy return node; } + isax_node *get_or_create_pRecBuf_root(parallel_first_buffer_layer *fbl, + root_mask_type mask, + isax_index *index, + pthread_mutex_t *lock_firstnode, + int total_workernumber) + { + parallel_fbl_soft_buffer *current_buffer = &fbl->soft_buffers[(int)mask]; + if (current_buffer->initialized) + return current_buffer->node; + + pthread_mutex_lock(lock_firstnode); + if (!current_buffer->initialized) + { + current_buffer->max_buffer_size = + (int *)calloc((size_t)total_workernumber, sizeof(int)); + current_buffer->buffer_size = + (int *)calloc((size_t)total_workernumber, sizeof(int)); + current_buffer->sax_records = + (sax_type **)calloc((size_t)total_workernumber, sizeof(sax_type *)); + current_buffer->pos_records = + (file_position_type **)calloc((size_t)total_workernumber, + sizeof(file_position_type *)); + current_buffer->node = + isax_root_node_init(mask, index->settings->initial_leaf_buffer_size); + + if (current_buffer->max_buffer_size == NULL || + current_buffer->buffer_size == NULL || + current_buffer->sax_records == NULL || + current_buffer->pos_records == NULL || + current_buffer->node == NULL) + { + free(current_buffer->max_buffer_size); + free(current_buffer->buffer_size); + free(current_buffer->sax_records); + free(current_buffer->pos_records); + if (current_buffer->node != NULL) + { + destroy_node_buffer(current_buffer->node->buffer); + free(current_buffer->node); + } + current_buffer->max_buffer_size = NULL; + current_buffer->buffer_size = NULL; + current_buffer->sax_records = NULL; + current_buffer->pos_records = NULL; + current_buffer->node = NULL; + pthread_mutex_unlock(lock_firstnode); + return NULL; + } + + current_buffer->node->is_leaf = 1; + current_buffer->node->previous = NULL; + current_buffer->node->next = index->first_node; + if (index->first_node != NULL) + index->first_node->previous = current_buffer->node; + index->first_node = current_buffer->node; + current_buffer->initialized = 1; + __sync_fetch_and_add(&(index->root_nodes), 1); + } + pthread_mutex_unlock(lock_firstnode); + return current_buffer->node; + } + isax_node *insert_to_pRecBuf(parallel_first_buffer_layer *fbl, sax_type *sax, file_position_type *pos, root_mask_type mask, isax_index *index, pthread_mutex_t *lock_firstnode, int workernumber, int total_workernumber) @@ -583,50 +649,10 @@ namespace daisy int current_buffer_number; // char *cd_s, *cd_p; - // Check if this buffer is initialized - - if (!current_buffer->initialized) - { - pthread_mutex_lock(lock_firstnode); - if (!current_buffer->initialized) - { - - current_buffer->max_buffer_size = (int *)malloc(sizeof(int) * total_workernumber); - current_buffer->buffer_size = (int *)malloc(sizeof(int) * total_workernumber); - current_buffer->sax_records = (sax_type **)malloc(sizeof(sax_type *) * total_workernumber); - current_buffer->pos_records = (file_position_type **)malloc(sizeof(file_position_type *) * total_workernumber); - for (int i = 0; i < total_workernumber; i++) - { - current_buffer->max_buffer_size[i] = 0; - current_buffer->buffer_size[i] = 0; - current_buffer->pos_records[i] = NULL; - current_buffer->sax_records[i] = NULL; - } - current_buffer->node = isax_root_node_init(mask, index->settings->initial_leaf_buffer_size); - current_buffer->node->is_leaf = 1; - current_buffer->initialized = 1; - if (index->first_node == NULL) - { - index->first_node = current_buffer->node; - pthread_mutex_unlock(lock_firstnode); - current_buffer->node->next = NULL; - current_buffer->node->previous = NULL; - } - else - { - isax_node *prev_first = index->first_node; - index->first_node = current_buffer->node; - index->first_node->next = prev_first; - prev_first->previous = current_buffer->node; - pthread_mutex_unlock(lock_firstnode); - } - __sync_fetch_and_add(&(index->root_nodes), 1); - } - else - { - pthread_mutex_unlock(lock_firstnode); - } - } + // Check if this buffer is initialized. + if (get_or_create_pRecBuf_root(fbl, mask, index, lock_firstnode, + total_workernumber) == NULL) + return NULL; // Check if this buffer is not full! if (current_buffer->buffer_size[workernumber] >= current_buffer->max_buffer_size[workernumber]) diff --git a/lib/isax/iSAXIndex.hpp b/lib/isax/iSAXIndex.hpp index e322fb8..65e4d09 100644 --- a/lib/isax/iSAXIndex.hpp +++ b/lib/isax/iSAXIndex.hpp @@ -407,6 +407,15 @@ namespace daisy isax_node *isax_leaf_node_init(int initial_buffer_size); isax_node *isax_root_node_init(root_mask_type mask, int initial_buffer_size); + // Return the root node for a SAX root mask, creating and registering it in + // the parallel first-buffer layer when it does not exist yet. This does not + // append a record to the FBL, so callers with stable record storage can + // insert directly into the tree without invalidating existing FBL pointers. + isax_node *get_or_create_pRecBuf_root(parallel_first_buffer_layer *fbl, + root_mask_type mask, + isax_index *index, + pthread_mutex_t *lock_firstnode, + int total_workernumber); isax_node *insert_to_pRecBuf(parallel_first_buffer_layer *fbl, sax_type *sax, file_position_type *pos, root_mask_type mask, isax_index *index, pthread_mutex_t *lock_firstnode, int workernumber, int total_workernumber); // EKOSMAS-specific versions for Odyssey @@ -485,4 +494,4 @@ namespace daisy } -#endif \ No newline at end of file +#endif diff --git a/pybinds/setup.cpp b/pybinds/setup.cpp index c781f0c..084820f 100644 --- a/pybinds/setup.cpp +++ b/pybinds/setup.cpp @@ -408,8 +408,12 @@ PYBIND11_MODULE(_core, m) .def("setWarpingWindow", &daisy::Messi::setWarpingWindow, "Set the warping window size for DTW") // Build the index from a 2D NumPy array - .def("buildIndex", [](daisy::Messi &self, pybind11::array_t db) + .def("buildIndex", [](daisy::Messi &self, pybind11::array db) { + if (!db.dtype().is(pybind11::dtype::of())) + throw std::runtime_error("Database array must have dtype float32"); + if ((db.flags() & pybind11::array::c_style) == 0) + throw std::runtime_error("Database array must be C-contiguous"); pybind11::buffer_info buf = db.request(); if (buf.ndim != 2) throw std::runtime_error("Database array must be 2D"); @@ -419,7 +423,35 @@ PYBIND11_MODULE(_core, m) // Create InMemoryDataSource from numpy array daisy::InMemoryDataSource data_source(static_cast(buf.ptr), n, d); - self.buildIndex(&data_source); }, "Build the MESSI index from a 2D float32 NumPy array") + self.buildIndex(&data_source); }, + pybind11::keep_alive<1, 2>(), + "Build the MESSI index from a contiguous 2D float32 NumPy array") + + // Streaming: append one series or a contiguous batch to the live tree. + .def("insert", [](daisy::Messi &self, + pybind11::array_t series) + { + pybind11::buffer_info buf = series.request(); + if (buf.ndim != 1) + throw std::runtime_error("insert expects a 1D float32 array"); + if (self.getDim() != 0 && static_cast(buf.shape[0]) != self.getDim()) + throw std::runtime_error("insert series dimension does not match the index dimension"); + self.insert(static_cast(buf.ptr)); }, + "Incrementally insert one series into the live MESSI index") + + .def("insertBatch", [](daisy::Messi &self, + pybind11::array_t batch) + { + pybind11::buffer_info buf = batch.request(); + if (buf.ndim != 2) + throw std::runtime_error("insertBatch expects a 2D float32 array"); + if (self.getDim() != 0 && static_cast(buf.shape[1]) != self.getDim()) + throw std::runtime_error("insertBatch series dimension does not match the index dimension"); + self.insertBatch(static_cast(buf.ptr), + static_cast(buf.shape[0])); }, + "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) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f3c6f08..517f0d1 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -358,6 +358,30 @@ gtest_discover_tests(test_LbBruteforce_Streaming WORKING_DIRECTORY ${CMAKE_SOURC message(STATUS "Tests discovered for test_Messi_L2Square.") endif() +# ////// MESSI streaming ////// +add_executable( + test_Messi_Streaming + test_Messi_Streaming.cpp + test_utils.cpp +) +target_link_libraries( + test_Messi_Streaming + PRIVATE + GTest::gtest_main + dino_lib + commons_lib + stdc++fs +) +target_include_directories(test_Messi_Streaming + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons +) +gtest_discover_tests( + test_Messi_Streaming + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} +) + # ////// ODYSSEY L2Square (conditional on MPI) ////// # Odyssey tests must run under mpirun (e.g. mpirun -np 4); one CTest entry runs all GTest cases. if(DEBUG_MSG) diff --git a/tests/test_Messi_Streaming.cpp b/tests/test_Messi_Streaming.cpp new file mode 100644 index 0000000..e996592 --- /dev/null +++ b/tests/test_Messi_Streaming.cpp @@ -0,0 +1,214 @@ +#include + +#include "../lib/algos/Bruteforce.hpp" +#include "../lib/algos/Messi.hpp" + +#include +#include +#include + +namespace +{ + constexpr int DIM = 32; + + std::vector makeSeries(int n, int first_phase = 0) + { + std::vector data(static_cast(n) * DIM); + for (int i = 0; i < n; ++i) + { + float *series = data.data() + static_cast(i) * DIM; + const float phase = static_cast(first_phase + i) * 0.29f; + double mean = 0.0; + for (int j = 0; j < DIM; ++j) + { + series[j] = std::sin(0.17f * j + phase) + + 0.45f * std::cos(0.41f * j - 0.3f * phase); + mean += series[j]; + } + mean /= DIM; + + double variance = 0.0; + for (int j = 0; j < DIM; ++j) + variance += (series[j] - mean) * (series[j] - mean); + const double stddev = std::sqrt(variance / DIM); + for (int j = 0; j < DIM; ++j) + series[j] = static_cast((series[j] - mean) / stddev); + } + return data; + } + + daisy::MessiConfig streamingConfig() + { + daisy::MessiConfig config; + config.search_workers = 1; + config.index_workers = 1; + config.leaf_size = 4; + config.paa_segments = 16; + return config; + } + + void expectL2MatchesBruteforce(daisy::Messi &search, + const std::vector &data, + int n_database, + const std::vector &queries, + int n_query, + int k) + { + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.setNumThreads(1); + ground_truth.buildIndex(const_cast(data.data()), n_database, DIM); + + std::vector expected_indices(static_cast(n_query) * k); + std::vector expected_distances(static_cast(n_query) * k); + std::vector actual_indices(static_cast(n_query) * k); + std::vector actual_distances(static_cast(n_query) * k); + + ground_truth.searchIndex(queries.data(), n_query, k, + expected_indices.data(), expected_distances.data()); + search.searchIndex(queries.data(), n_query, k, + actual_indices.data(), actual_distances.data()); + + for (size_t i = 0; i < actual_indices.size(); ++i) + { + EXPECT_EQ(actual_indices[i], expected_indices[i]); + EXPECT_NEAR(actual_distances[i], expected_distances[i], 1e-4f); + } + } +} + +TEST(MessiStreamingTest, RequiresBuildAndRejectsNullInput) +{ + daisy::Messi search(daisy::DistanceType::L2_SQUARED, streamingConfig()); + EXPECT_THROW(search.insert(nullptr), std::runtime_error); + + auto initial = makeSeries(2); + search.buildIndex(initial.data(), 2, DIM); + EXPECT_THROW(search.insert(nullptr), std::invalid_argument); + EXPECT_NO_THROW(search.insertBatch(nullptr, 0)); +} + +TEST(MessiStreamingTest, SingleAndBatchInsertsMatchBruteforce) +{ + auto all = makeSeries(18); + auto queries = makeSeries(5, 50); + daisy::Messi search(daisy::DistanceType::L2_SQUARED, streamingConfig()); + search.buildIndex(all.data(), 7, DIM); + expectL2MatchesBruteforce(search, all, 7, queries, 5, 4); + + daisy::SimilaritySearchAlgorithm *streaming = &search; + streaming->insert(all.data() + 7 * DIM); + expectL2MatchesBruteforce(search, all, 8, queries, 5, 4); + + streaming->insertBatch(all.data() + 8 * DIM, 10); + ASSERT_EQ(search.getNDatabase(), 18u); + expectL2MatchesBruteforce(search, all, 18, queries, 5, 4); + + daisy::idx_t index = 0; + float distance = -1.0f; + search.searchIndex(all.data() + 17 * DIM, 1, 1, &index, &distance); + EXPECT_EQ(index, 17u); + EXPECT_FLOAT_EQ(distance, 0.0f); + + daisy::SearchConfig range; + range.type = daisy::QueryType::RANGE; + range.r = 1e-5f; + std::vector> indices; + std::vector> distances; + search.searchIndex(all.data() + 7 * DIM, 1, range, indices, distances); + ASSERT_EQ(indices.size(), 1u); + ASSERT_EQ(indices[0].size(), 1u); + EXPECT_EQ(indices[0][0], 7u); + EXPECT_FLOAT_EQ(distances[0][0], 0.0f); +} + +TEST(MessiStreamingTest, CreatesMissingRootAndSplitsExistingLeaf) +{ + std::vector initial(static_cast(2) * DIM); + for (int row = 0; row < 2; ++row) + { + for (int j = 0; j < DIM; ++j) + { + const int segment = j / 8; + initial[static_cast(row) * DIM + j] = + static_cast(segment * 2 - 3); + } + } + std::vector opposite(DIM); + for (int j = 0; j < DIM; ++j) + opposite[j] = -initial[j]; + + daisy::MessiConfig config; + config.search_workers = 1; + config.index_workers = 1; + config.leaf_size = 2; + config.paa_segments = 4; + daisy::Messi search(daisy::DistanceType::L2_SQUARED, config); + search.buildIndex(initial.data(), 2, DIM); + + ASSERT_NE(search.getIndex(), nullptr); + ASSERT_NE(search.getIndex()->first_node, nullptr); + daisy::isax_node *initial_root = search.getIndex()->first_node; + const unsigned long roots_before = search.getIndex()->root_nodes; + + // Same SAX root at capacity: the incremental insert must split the leaf. + search.insert(initial.data()); + EXPECT_FALSE(initial_root->is_leaf); + + // Negating every segment flips the root SAX mask, forcing root creation. + search.insert(opposite.data()); + EXPECT_GT(search.getIndex()->root_nodes, roots_before); + + daisy::idx_t index = 0; + float distance = -1.0f; + search.searchIndex(opposite.data(), 1, 1, &index, &distance); + EXPECT_EQ(index, 3u); + EXPECT_FLOAT_EQ(distance, 0.0f); +} + +TEST(MessiStreamingTest, DtwSearchIncludesInsertedSeries) +{ + auto all = makeSeries(8, 80); + daisy::Messi search(daisy::DistanceType::DTW, streamingConfig()); + search.buildIndex(all.data(), 4, DIM); + search.insert(all.data() + 4 * DIM); + search.insertBatch(all.data() + 5 * DIM, 3); + + daisy::idx_t index = 0; + float distance = -1.0f; + search.searchIndex(all.data() + 7 * DIM, 1, 1, &index, &distance); + EXPECT_EQ(index, 7u); + EXPECT_NEAR(distance, 0.0f, 1e-6f); +} + +TEST(MessiStreamingTest, EquidepthInsertsReuseInitialBreakpoints) +{ + auto all = makeSeries(16, 110); + auto queries = makeSeries(4, 140); + for (float &value : all) + value = 30.0f + 6.0f * value; + for (float &value : queries) + value = 30.0f + 6.0f * value; + + daisy::Messi search(daisy::DistanceType::L2_SQUARED, streamingConfig()); + search.setNormalized(false); + search.buildIndex(all.data(), 8, DIM); + search.insertBatch(all.data() + 8 * DIM, 8); + + expectL2MatchesBruteforce(search, all, 16, queries, 4, 4); +} + +TEST(MessiStreamingTest, CanInsertFromItsOwnDatabaseAcrossReallocation) +{ + auto initial = makeSeries(4, 180); + daisy::Messi search(daisy::DistanceType::L2_SQUARED, streamingConfig()); + search.buildIndex(initial.data(), 4, DIM); + + search.insert(search.getDatabase() + DIM); + ASSERT_EQ(search.getNDatabase(), 5u); + const float *owned_database = search.getDatabase(); + search.insertBatch(owned_database, 4); + ASSERT_EQ(search.getNDatabase(), 9u); + + for (int i = 0; i < 4 * DIM; ++i) + EXPECT_FLOAT_EQ(search.getDatabase()[5 * DIM + i], initial[i]); +} diff --git a/tests/test_streaming.py b/tests/test_streaming.py index ae0ea65..2035146 100644 --- a/tests/test_streaming.py +++ b/tests/test_streaming.py @@ -1,18 +1,23 @@ +import gc import unittest +import weakref import numpy as np -from daisy import BruteForceSearch, DistanceType, LbBruteforce +from daisy import BruteForceSearch, DistanceType, LbBruteforce, Messi class StreamingBindingsTest(unittest.TestCase): - def test_bruteforce_and_lb_bruteforce(self): + def test_streaming_algorithms(self): rng = np.random.default_rng(123) data = rng.normal(size=(8, 32)).astype(np.float32) - for algorithm in (BruteForceSearch, LbBruteforce): + for algorithm in (BruteForceSearch, LbBruteforce, Messi): with self.subTest(algorithm=algorithm.__name__): index = algorithm(DistanceType.L2_SQUARED) + if algorithm is Messi: + index.setIndexWorkers(1) + index.setSearchWorkers(1) index.buildIndex(data[:3]) index.insert(data[3]) index.insertBatch(data[4:]) @@ -22,15 +27,54 @@ def test_bruteforce_and_lb_bruteforce(self): self.assertAlmostEqual(float(distances[0, 0]), 0.0, places=6) def test_dimension_validation(self): - index = BruteForceSearch(DistanceType.L2_SQUARED) - index.buildIndex(np.zeros((2, 32), dtype=np.float32)) + for algorithm in (BruteForceSearch, LbBruteforce, Messi): + with self.subTest(algorithm=algorithm.__name__): + index = algorithm(DistanceType.L2_SQUARED) + if algorithm is Messi: + index.setIndexWorkers(1) + index.setSearchWorkers(1) + index.buildIndex(np.zeros((2, 32), dtype=np.float32)) + + with self.assertRaises(RuntimeError): + index.insert(np.zeros(31, dtype=np.float32)) + with self.assertRaises(RuntimeError): + index.insertBatch(np.zeros((2, 31), dtype=np.float32)) + + def test_messi_keeps_borrowed_build_array_alive(self): + rng = np.random.default_rng(456) + + def build_from_temporary(): + data = rng.normal(size=(5, 32)).astype(np.float32) + index = Messi(DistanceType.L2_SQUARED) + index.setIndexWorkers(1) + index.setSearchWorkers(1) + borrowed = data[:3].copy() + borrowed_ref = weakref.ref(borrowed) + index.buildIndex(borrowed) + return index, borrowed_ref, data[0].copy(), data[3].copy() + + index, borrowed_ref, initial, inserted = build_from_temporary() + gc.collect() + self.assertIsNotNone(borrowed_ref()) + initial_indices, initial_distances = index.searchIndex(initial.reshape(1, -1), 1) + self.assertEqual(int(initial_indices[0, 0]), 0) + self.assertAlmostEqual(float(initial_distances[0, 0]), 0.0, places=6) + + index.insert(inserted) + indices, distances = index.searchIndex(inserted.reshape(1, -1), 1) + self.assertEqual(int(indices[0, 0]), 3) + self.assertAlmostEqual(float(distances[0, 0]), 0.0, places=6) + + def test_messi_build_rejects_temporary_numpy_conversions(self): + index = Messi(DistanceType.L2_SQUARED) + index.setIndexWorkers(1) + index.setSearchWorkers(1) with self.assertRaises(RuntimeError): - index.insert(np.zeros(31, dtype=np.float32)) + index.buildIndex(np.zeros((3, 32), dtype=np.float64)) with self.assertRaises(RuntimeError): - index.insertBatch(np.zeros((2, 31), dtype=np.float32)) + index.buildIndex(np.zeros((3, 64), dtype=np.float32)[:, ::2]) if __name__ == "__main__": unittest.main() - From 72a0d6601c6d0d9e49b50e51eb866eaa0972b0b1 Mon Sep 17 00:00:00 2001 From: sophisid Date: Thu, 17 Sep 2026 11:22:02 +0200 Subject: [PATCH 2/4] fix minor issues --- demos/demo_Messi_Streaming.cpp | 27 +++++++++++++------- demos/demo_Messi_Streaming.py | 19 ++++++++++---- lib/algos/Messi.cpp | 10 ++++++-- tests/test_Messi_Streaming.cpp | 46 ++++++++++++++++++++++++++++++---- tests/test_streaming.py | 11 +++++++- 5 files changed, 91 insertions(+), 22 deletions(-) diff --git a/demos/demo_Messi_Streaming.cpp b/demos/demo_Messi_Streaming.cpp index 40fc336..1822c70 100644 --- a/demos/demo_Messi_Streaming.cpp +++ b/demos/demo_Messi_Streaming.cpp @@ -21,18 +21,27 @@ int main() search.setSearchWorkers(4); search.buildIndex(stream, initial, dim); - search.insert(stream + initial * dim); - search.insertBatch(stream + (initial + 1) * dim, batch - 1); - daisy::idx_t *indices = new daisy::idx_t[n_query * k]; float *distances = new float[n_query * k]; - search.searchIndex(query, n_query, k, indices, distances); - std::printf("MESSI now contains %llu series. Query 0 kNN: ", - search.getNDatabase()); - for (daisy::idx_t j = 0; j < k; ++j) - std::printf("%llu(%.3f) ", indices[j], distances[j]); - std::printf("\n"); + // The index stays queryable between updates: search after every addition. + auto queryAndReport = [&](const char *stage) + { + search.searchIndex(query, n_query, k, indices, distances); + std::printf("%-16s MESSI contains %llu series. Query 0 kNN: ", + stage, search.getNDatabase()); + for (daisy::idx_t j = 0; j < k; ++j) + std::printf("%llu(%.3f) ", indices[j], distances[j]); + std::printf("\n"); + }; + + queryAndReport("after build"); + + search.insert(stream + initial * dim); + queryAndReport("after insert"); + + search.insertBatch(stream + (initial + 1) * dim, batch - 1); + queryAndReport("after batch"); delete[] stream; delete[] query; diff --git a/demos/demo_Messi_Streaming.py b/demos/demo_Messi_Streaming.py index e416160..7d1b2d9 100644 --- a/demos/demo_Messi_Streaming.py +++ b/demos/demo_Messi_Streaming.py @@ -14,10 +14,19 @@ index.setSearchWorkers(4) index.buildIndex(stream[:5000]) + +def query_and_report(stage, size): + """The index stays queryable between updates: search after every addition.""" + indices, distances = index.searchIndex(queries, 5) + print(f"{stage}: MESSI size {size}") + print(" Query 0 IDs:", indices[0]) + print(" Query 0 distances:", distances[0]) + + +query_and_report("after build", 5000) + index.insert(stream[5000]) -index.insertBatch(stream[5001:]) +query_and_report("after insert", 5001) -indices, distances = index.searchIndex(queries, 5) -print("MESSI size:", 6000) -print("Query 0 IDs:", indices[0]) -print("Query 0 distances:", distances[0]) +index.insertBatch(stream[5001:]) +query_and_report("after batch", 6000) diff --git a/lib/algos/Messi.cpp b/lib/algos/Messi.cpp index 345ba68..4134c71 100644 --- a/lib/algos/Messi.cpp +++ b/lib/algos/Messi.cpp @@ -1088,6 +1088,11 @@ namespace daisy void Messi::searchIndexL2Squared(const float *query, const idx_t n_query, const idx_t k, idx_t *I, float *D) { + // Lower bounds read the active breakpoints, which are global: interleaving + // searches with inserts (or with another live index) can leave a different + // table installed, so reinstall ours before every search. + activateBreakpoints(); + ts_type *paa = (ts_type *)malloc(sizeof(ts_type) * index->settings->paa_segments); node_list nodelist; @@ -1162,6 +1167,8 @@ namespace daisy void Messi::searchIndexDTW(const float *query, const idx_t n_query, const idx_t k, idx_t *I, float *D) { + activateBreakpoints(); + isax_index *index = this->index; node_list nodelist; @@ -1295,13 +1302,13 @@ namespace daisy std::vector> &I, std::vector> &D) { - activateBreakpoints(); if (config.type == QueryType::TOP_K) { SimilaritySearchAlgorithm::searchIndex(query, n_query, config, I, D); return; } + activateBreakpoints(); ts_type *paa = (ts_type *)malloc(sizeof(ts_type) * index->settings->paa_segments); node_list nodelist; nodelist.nlist = (isax_node **)malloc(sizeof(isax_node *) * (int)pow(2, index->settings->paa_segments)); @@ -1339,7 +1346,6 @@ namespace daisy void Messi::searchIndex(const float *query, const idx_t n_query, const idx_t k, idx_t *I, float *D) { - activateBreakpoints(); if (this->distance_type == DistanceType::L2_SQUARED) { searchIndexL2Squared(query, n_query, k, I, D); diff --git a/tests/test_Messi_Streaming.cpp b/tests/test_Messi_Streaming.cpp index e996592..9bbc220 100644 --- a/tests/test_Messi_Streaming.cpp +++ b/tests/test_Messi_Streaming.cpp @@ -150,16 +150,22 @@ TEST(MessiStreamingTest, CreatesMissingRootAndSplitsExistingLeaf) daisy::isax_node *initial_root = search.getIndex()->first_node; const unsigned long roots_before = search.getIndex()->root_nodes; + daisy::idx_t index = 0; + float distance = -1.0f; + // Same SAX root at capacity: the incremental insert must split the leaf. search.insert(initial.data()); EXPECT_FALSE(initial_root->is_leaf); + // The split tree stays searchable before the next insert arrives. + search.searchIndex(initial.data(), 1, 1, &index, &distance); + EXPECT_LT(index, search.getNDatabase()); + EXPECT_FLOAT_EQ(distance, 0.0f); + // Negating every segment flips the root SAX mask, forcing root creation. search.insert(opposite.data()); EXPECT_GT(search.getIndex()->root_nodes, roots_before); - daisy::idx_t index = 0; - float distance = -1.0f; search.searchIndex(opposite.data(), 1, 1, &index, &distance); EXPECT_EQ(index, 3u); EXPECT_FLOAT_EQ(distance, 0.0f); @@ -170,11 +176,19 @@ TEST(MessiStreamingTest, DtwSearchIncludesInsertedSeries) auto all = makeSeries(8, 80); daisy::Messi search(daisy::DistanceType::DTW, streamingConfig()); search.buildIndex(all.data(), 4, DIM); - search.insert(all.data() + 4 * DIM); - search.insertBatch(all.data() + 5 * DIM, 3); daisy::idx_t index = 0; float distance = -1.0f; + search.searchIndex(all.data() + 3 * DIM, 1, 1, &index, &distance); + EXPECT_EQ(index, 3u); + EXPECT_NEAR(distance, 0.0f, 1e-6f); + + search.insert(all.data() + 4 * DIM); + search.searchIndex(all.data() + 4 * DIM, 1, 1, &index, &distance); + EXPECT_EQ(index, 4u); + EXPECT_NEAR(distance, 0.0f, 1e-6f); + + search.insertBatch(all.data() + 5 * DIM, 3); search.searchIndex(all.data() + 7 * DIM, 1, 1, &index, &distance); EXPECT_EQ(index, 7u); EXPECT_NEAR(distance, 0.0f, 1e-6f); @@ -192,8 +206,12 @@ TEST(MessiStreamingTest, EquidepthInsertsReuseInitialBreakpoints) daisy::Messi search(daisy::DistanceType::L2_SQUARED, streamingConfig()); search.setNormalized(false); search.buildIndex(all.data(), 8, DIM); - search.insertBatch(all.data() + 8 * DIM, 8); + expectL2MatchesBruteforce(search, all, 8, queries, 4, 4); + + search.insertBatch(all.data() + 8 * DIM, 4); + expectL2MatchesBruteforce(search, all, 12, queries, 4, 4); + search.insertBatch(all.data() + 12 * DIM, 4); expectL2MatchesBruteforce(search, all, 16, queries, 4, 4); } @@ -205,10 +223,28 @@ TEST(MessiStreamingTest, CanInsertFromItsOwnDatabaseAcrossReallocation) search.insert(search.getDatabase() + DIM); ASSERT_EQ(search.getNDatabase(), 5u); + + // The series now appears twice, so assert on the hit contents instead of its id. + daisy::idx_t index = 0; + float distance = -1.0f; + search.searchIndex(initial.data() + DIM, 1, 1, &index, &distance); + EXPECT_FLOAT_EQ(distance, 0.0f); + ASSERT_LT(index, search.getNDatabase()); + for (int j = 0; j < DIM; ++j) + EXPECT_FLOAT_EQ(search.getDatabase()[static_cast(index) * DIM + j], + initial[DIM + j]); + const float *owned_database = search.getDatabase(); search.insertBatch(owned_database, 4); ASSERT_EQ(search.getNDatabase(), 9u); for (int i = 0; i < 4 * DIM; ++i) EXPECT_FLOAT_EQ(search.getDatabase()[5 * DIM + i], initial[i]); + + search.searchIndex(initial.data() + 3 * DIM, 1, 1, &index, &distance); + EXPECT_FLOAT_EQ(distance, 0.0f); + ASSERT_LT(index, search.getNDatabase()); + for (int j = 0; j < DIM; ++j) + EXPECT_FLOAT_EQ(search.getDatabase()[static_cast(index) * DIM + j], + initial[3 * DIM + j]); } diff --git a/tests/test_streaming.py b/tests/test_streaming.py index 2035146..c1f79bf 100644 --- a/tests/test_streaming.py +++ b/tests/test_streaming.py @@ -19,9 +19,18 @@ def test_streaming_algorithms(self): index.setIndexWorkers(1) index.setSearchWorkers(1) index.buildIndex(data[:3]) + + # Every addition leaves the index queryable, so search in between. + indices, distances = index.searchIndex(data[2:3], 1) + self.assertEqual(int(indices[0, 0]), 2) + self.assertAlmostEqual(float(distances[0, 0]), 0.0, places=6) + index.insert(data[3]) - index.insertBatch(data[4:]) + indices, distances = index.searchIndex(data[3:4], 1) + self.assertEqual(int(indices[0, 0]), 3) + self.assertAlmostEqual(float(distances[0, 0]), 0.0, places=6) + index.insertBatch(data[4:]) indices, distances = index.searchIndex(data[7:8], 1) self.assertEqual(int(indices[0, 0]), 7) self.assertAlmostEqual(float(distances[0, 0]), 0.0, places=6) From ee23d2eab0e0914f528d39e279405c4b6fc697fa Mon Sep 17 00:00:00 2001 From: sophisid Date: Thu, 17 Sep 2026 11:47:41 +0200 Subject: [PATCH 3/4] range queries demos --- README.md | 19 ++++ demos/CMakeLists.txt | 77 +++++++++++++ demos/demo_Bruteforce_Range.cpp | 62 ++++++++++ demos/demo_Coconut_Range.cpp | 98 ++++++++++++++++ demos/demo_DumpyOS_Range.cpp | 98 ++++++++++++++++ demos/demo_Fresh_Range.cpp | 98 ++++++++++++++++ demos/demo_Hercules_Range.cpp | 102 +++++++++++++++++ demos/demo_LbBruteforce_Range.cpp | 98 ++++++++++++++++ demos/demo_Messi_Range.cpp | 98 ++++++++++++++++ demos/demo_Odyssey_Range.cpp | 182 ++++++++++++++++++++++++++++++ demos/demo_ParIS_Range.cpp | 112 ++++++++++++++++++ demos/demo_Sing_Range.cpp | 110 ++++++++++++++++++ demos/demo_Sofa_Range.cpp | 98 ++++++++++++++++ docs/demos-guide.md | 10 +- 14 files changed, 1261 insertions(+), 1 deletion(-) create mode 100644 demos/demo_Bruteforce_Range.cpp create mode 100644 demos/demo_Coconut_Range.cpp create mode 100644 demos/demo_DumpyOS_Range.cpp create mode 100644 demos/demo_Fresh_Range.cpp create mode 100644 demos/demo_Hercules_Range.cpp create mode 100644 demos/demo_LbBruteforce_Range.cpp create mode 100644 demos/demo_Messi_Range.cpp create mode 100644 demos/demo_Odyssey_Range.cpp create mode 100644 demos/demo_ParIS_Range.cpp create mode 100644 demos/demo_Sing_Range.cpp create mode 100644 demos/demo_Sofa_Range.cpp diff --git a/README.md b/README.md index b7e012d..f4754ec 100644 --- a/README.md +++ b/README.md @@ -79,6 +79,25 @@ index. Its first insert copies a borrowed initial in-memory database into owned Streaming updates are not concurrent with queries. Inserts can reallocate the owned database, so callers should not retain a pointer returned by `getDatabase()` across them. +### Range (distance-r) queries + +Every algorithm except Coconut's streaming-only paths answers range queries through +`SearchConfig`. Instead of a fixed `k`, each query returns however many series fall within the +radius, so the results come back as one vector per query: + +```cpp +daisy::SearchConfig config; +config.type = daisy::QueryType::RANGE; +config.r = radius; // squared L2 distance + +std::vector> I; +std::vector> D; +search.searchIndex(query, n_query, config, I, D); +``` + +`demos/demo__Range.cpp` shows this for each algorithm and cross-checks the returned +sets against brute force. + ## Quickstart diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt index 209426c..6e04dbf 100644 --- a/demos/CMakeLists.txt +++ b/demos/CMakeLists.txt @@ -44,6 +44,13 @@ if(BUILD_DEMO) ${CMAKE_CURRENT_SOURCE_DIR}/../commons ) + add_executable(demo_Bruteforce_Range demo_Bruteforce_Range.cpp) + target_link_libraries(demo_Bruteforce_Range PRIVATE dino_lib commons_lib) + target_include_directories(demo_Bruteforce_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + # ////// COCONUT (static + streaming) ////// if(BUILD_COCONUT) if(DEBUG_MSG) @@ -61,6 +68,13 @@ if(BUILD_DEMO) ${CMAKE_CURRENT_SOURCE_DIR}/../lib ${CMAKE_CURRENT_SOURCE_DIR}/../commons ) + + add_executable(demo_Coconut_Range demo_Coconut_Range.cpp) + target_link_libraries(demo_Coconut_Range PRIVATE dino_lib commons_lib) + target_include_directories(demo_Coconut_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) elseif(DEBUG_MSG) message(STATUS "BUILD_COCONUT is OFF. Skipping Coconut demos.") endif() @@ -138,6 +152,13 @@ if(BUILD_DEMO) ${CMAKE_CURRENT_SOURCE_DIR}/../commons ) + add_executable(demo_LbBruteforce_Range demo_LbBruteforce_Range.cpp) + target_link_libraries(demo_LbBruteforce_Range PRIVATE dino_lib commons_lib) + target_include_directories(demo_LbBruteforce_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + # ////// LBBRUTEFORCE DTW ////// if(DEBUG_MSG) message(STATUS "---") @@ -211,6 +232,13 @@ if(BUILD_DEMO) ${CMAKE_CURRENT_SOURCE_DIR}/../commons ) + add_executable(demo_Messi_Range demo_Messi_Range.cpp) + target_link_libraries(demo_Messi_Range PRIVATE dino_lib commons_lib) + target_include_directories(demo_Messi_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + # ////// MESSI DTW ////// if(DEBUG_MSG) message(STATUS "---") @@ -314,6 +342,13 @@ if(BUILD_DEMO) ${CMAKE_CURRENT_SOURCE_DIR}/../lib ${CMAKE_CURRENT_SOURCE_DIR}/../commons ) + + add_executable(demo_Odyssey_Range demo_Odyssey_Range.cpp) + target_link_libraries(demo_Odyssey_Range PRIVATE dino_lib commons_lib MPI::MPI_CXX) + target_include_directories(demo_Odyssey_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) else() if(DEBUG_MSG) message(STATUS "BUILD_ODYSSEY_AVAILABLE is FALSE. demo_Odyssey_L2Square (MPI-dependent) will NOT be built.") @@ -386,6 +421,13 @@ if(BUILD_DEMO) message(STATUS "Include directories added for demo_ParIS_DTW.") endif() + add_executable(demo_ParIS_Range demo_ParIS_Range.cpp) + target_link_libraries(demo_ParIS_Range PRIVATE dino_lib commons_lib) + target_include_directories(demo_ParIS_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + # ////// SING L2Square ////// if(DEBUG_MSG) message(STATUS "---") @@ -403,6 +445,13 @@ if(BUILD_DEMO) ${CMAKE_CURRENT_SOURCE_DIR}/../lib ${CMAKE_CURRENT_SOURCE_DIR}/../commons ) + + add_executable(demo_Sing_Range demo_Sing_Range.cpp) + target_link_libraries(demo_Sing_Range PRIVATE dino_lib commons_lib) + target_include_directories(demo_Sing_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) else() if(DEBUG_MSG) message(STATUS "BUILD_SING_AVAILABLE is FALSE. demo_Sing_L2Square will NOT be built.") @@ -442,6 +491,13 @@ if(BUILD_DEMO) if(DEBUG_MSG) message(STATUS "Include directories added for demo_Sofa_L2Square.") endif() + + add_executable(demo_Sofa_Range demo_Sofa_Range.cpp) + target_link_libraries(demo_Sofa_Range PRIVATE dino_lib commons_lib) + target_include_directories(demo_Sofa_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) else() if(DEBUG_MSG) message(STATUS "BUILD_SOFA_AVAILABLE is FALSE. demo_Sofa_L2Square (FFTW3-dependent) will NOT be built.") @@ -478,6 +534,13 @@ if(BUILD_DEMO) message(STATUS "Include directories added for demo_Hercules_L2Square.") endif() + add_executable(demo_Hercules_Range demo_Hercules_Range.cpp) + target_link_libraries(demo_Hercules_Range PRIVATE dino_lib commons_lib) + target_include_directories(demo_Hercules_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + # ////// DUMPYOS L2Square ////// if(DEBUG_MSG) message(STATUS "---") @@ -538,6 +601,13 @@ if(BUILD_DEMO) message(STATUS "Include directories added for demo_DumpyOS_DTW.") endif() + add_executable(demo_DumpyOS_Range demo_DumpyOS_Range.cpp) + target_link_libraries(demo_DumpyOS_Range PRIVATE dino_lib commons_lib) + target_include_directories(demo_DumpyOS_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + # ////// FRESH L2Square ////// if(DEBUG_MSG) message(STATUS "---") @@ -598,6 +668,13 @@ if(BUILD_DEMO) message(STATUS "Include directories added for demo_Fresh_DTW.") endif() + add_executable(demo_Fresh_Range demo_Fresh_Range.cpp) + target_link_libraries(demo_Fresh_Range PRIVATE dino_lib commons_lib) + target_include_directories(demo_Fresh_Range PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + else() if(DEBUG_MSG) diff --git a/demos/demo_Bruteforce_Range.cpp b/demos/demo_Bruteforce_Range.cpp new file mode 100644 index 0000000..68913ed --- /dev/null +++ b/demos/demo_Bruteforce_Range.cpp @@ -0,0 +1,62 @@ +// Bruteforce range search: report every series within squared-L2 distance r of a query. +// Brute force is the exact baseline that every other range demo checks itself against. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + float *database = loadRandomData(n_database, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + printf("Loaded %llu database points and %llu query points with dimension %llu\n", + n_database, n_query, dim); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + daisy::BruteForceSearch bf_search(daisy::DistanceType::L2_SQUARED); + bf_search.setNumThreads(4); + bf_search.buildIndex(database, n_database, dim); + + std::vector> I; + std::vector> D; + bf_search.searchIndex(query, n_query, config, I, D); + reportRange("Bruteforce", r, I, D); + + delete[] database; + delete[] query; + + return 0; +} diff --git a/demos/demo_Coconut_Range.cpp b/demos/demo_Coconut_Range.cpp new file mode 100644 index 0000000..29dda7b --- /dev/null +++ b/demos/demo_Coconut_Range.cpp @@ -0,0 +1,98 @@ +// Coconut range search: report every series within squared-L2 distance r of a query. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +// Range search returns an unordered set of hits, so compare sets rather than positions. +// "missing" are hits brute force found and the index did not; "extra" are the other way +// round, which for an exact index means false positives. +static bool matchesBruteforce(const char *name, + const std::vector> &I, + const std::vector> >_I) +{ + bool all_match = true; + printf("\n=== Verification against brute force ===\n"); + for (size_t qi = 0; qi < I.size(); qi++) + { + const std::set got(I[qi].begin(), I[qi].end()); + const std::set expected(gt_I[qi].begin(), gt_I[qi].end()); + + size_t missing = 0; + for (daisy::idx_t idx : expected) + if (got.count(idx) == 0) missing++; + size_t extra = 0; + for (daisy::idx_t idx : got) + if (expected.count(idx) == 0) extra++; + + all_match = all_match && (missing == 0 && extra == 0); + printf("Query %zu: %s=%zu, bruteforce=%zu, missing=%zu, extra=%zu\n", + qi, name, got.size(), expected.size(), missing, extra); + } + return all_match; +} + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + float *database = loadRandomData(n_database, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + printf("Loaded %llu database points and %llu query points with dimension %llu\n", + n_database, n_query, dim); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + daisy::Coconut coconut(daisy::DistanceType::L2_SQUARED); + coconut.setNumThreads(4); + coconut.buildIndex(database, n_database, dim); + + std::vector> I; + std::vector> D; + coconut.searchIndex(query, n_query, config, I, D); + reportRange("Coconut", r, I, D); + + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.buildIndex(database, n_database, dim); + std::vector> gt_I; + std::vector> gt_D; + ground_truth.searchIndex(query, n_query, config, gt_I, gt_D); + + printf("\nAll queries match: %s\n", + matchesBruteforce("Coconut", I, gt_I) ? "yes" : "no"); + + delete[] database; + delete[] query; + + return 0; +} diff --git a/demos/demo_DumpyOS_Range.cpp b/demos/demo_DumpyOS_Range.cpp new file mode 100644 index 0000000..c4d8e64 --- /dev/null +++ b/demos/demo_DumpyOS_Range.cpp @@ -0,0 +1,98 @@ +// DumpyOS range search: report every series within squared-L2 distance r of a query. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +// Range search returns an unordered set of hits, so compare sets rather than positions. +// "missing" are hits brute force found and the index did not; "extra" are the other way +// round, which for an exact index means false positives. +static bool matchesBruteforce(const char *name, + const std::vector> &I, + const std::vector> >_I) +{ + bool all_match = true; + printf("\n=== Verification against brute force ===\n"); + for (size_t qi = 0; qi < I.size(); qi++) + { + const std::set got(I[qi].begin(), I[qi].end()); + const std::set expected(gt_I[qi].begin(), gt_I[qi].end()); + + size_t missing = 0; + for (daisy::idx_t idx : expected) + if (got.count(idx) == 0) missing++; + size_t extra = 0; + for (daisy::idx_t idx : got) + if (expected.count(idx) == 0) extra++; + + all_match = all_match && (missing == 0 && extra == 0); + printf("Query %zu: %s=%zu, bruteforce=%zu, missing=%zu, extra=%zu\n", + qi, name, got.size(), expected.size(), missing, extra); + } + return all_match; +} + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + float *database = loadRandomData(n_database, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + printf("Loaded %llu database points and %llu query points with dimension %llu\n", + n_database, n_query, dim); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + daisy::DumpyOS dumpyos_search(daisy::DistanceType::L2_SQUARED); + dumpyos_search.setNumThreads(4); + dumpyos_search.buildIndex(database, n_database, dim); + + std::vector> I; + std::vector> D; + dumpyos_search.searchIndex(query, n_query, config, I, D); + reportRange("DumpyOS", r, I, D); + + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.buildIndex(database, n_database, dim); + std::vector> gt_I; + std::vector> gt_D; + ground_truth.searchIndex(query, n_query, config, gt_I, gt_D); + + printf("\nAll queries match: %s\n", + matchesBruteforce("DumpyOS", I, gt_I) ? "yes" : "no"); + + delete[] database; + delete[] query; + + return 0; +} diff --git a/demos/demo_Fresh_Range.cpp b/demos/demo_Fresh_Range.cpp new file mode 100644 index 0000000..d026cc7 --- /dev/null +++ b/demos/demo_Fresh_Range.cpp @@ -0,0 +1,98 @@ +// Fresh range search: report every series within squared-L2 distance r of a query. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +// Range search returns an unordered set of hits, so compare sets rather than positions. +// "missing" are hits brute force found and the index did not; "extra" are the other way +// round, which for an exact index means false positives. +static bool matchesBruteforce(const char *name, + const std::vector> &I, + const std::vector> >_I) +{ + bool all_match = true; + printf("\n=== Verification against brute force ===\n"); + for (size_t qi = 0; qi < I.size(); qi++) + { + const std::set got(I[qi].begin(), I[qi].end()); + const std::set expected(gt_I[qi].begin(), gt_I[qi].end()); + + size_t missing = 0; + for (daisy::idx_t idx : expected) + if (got.count(idx) == 0) missing++; + size_t extra = 0; + for (daisy::idx_t idx : got) + if (expected.count(idx) == 0) extra++; + + all_match = all_match && (missing == 0 && extra == 0); + printf("Query %zu: %s=%zu, bruteforce=%zu, missing=%zu, extra=%zu\n", + qi, name, got.size(), expected.size(), missing, extra); + } + return all_match; +} + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + float *database = loadRandomData(n_database, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + printf("Loaded %llu database points and %llu query points with dimension %llu\n", + n_database, n_query, dim); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + daisy::Fresh fresh_search(daisy::DistanceType::L2_SQUARED); + fresh_search.setNumThreads(4); + fresh_search.buildIndex(database, n_database, dim); + + std::vector> I; + std::vector> D; + fresh_search.searchIndex(query, n_query, config, I, D); + reportRange("Fresh", r, I, D); + + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.buildIndex(database, n_database, dim); + std::vector> gt_I; + std::vector> gt_D; + ground_truth.searchIndex(query, n_query, config, gt_I, gt_D); + + printf("\nAll queries match: %s\n", + matchesBruteforce("Fresh", I, gt_I) ? "yes" : "no"); + + delete[] database; + delete[] query; + + return 0; +} diff --git a/demos/demo_Hercules_Range.cpp b/demos/demo_Hercules_Range.cpp new file mode 100644 index 0000000..bc4073d --- /dev/null +++ b/demos/demo_Hercules_Range.cpp @@ -0,0 +1,102 @@ +// Hercules range search: report every series within squared-L2 distance r of a query. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +// Range search returns an unordered set of hits, so compare sets rather than positions. +// "missing" are hits brute force found and the index did not; "extra" are the other way +// round, which for an exact index means false positives. +static bool matchesBruteforce(const char *name, + const std::vector> &I, + const std::vector> >_I) +{ + bool all_match = true; + printf("\n=== Verification against brute force ===\n"); + for (size_t qi = 0; qi < I.size(); qi++) + { + const std::set got(I[qi].begin(), I[qi].end()); + const std::set expected(gt_I[qi].begin(), gt_I[qi].end()); + + size_t missing = 0; + for (daisy::idx_t idx : expected) + if (got.count(idx) == 0) missing++; + size_t extra = 0; + for (daisy::idx_t idx : got) + if (expected.count(idx) == 0) extra++; + + all_match = all_match && (missing == 0 && extra == 0); + printf("Query %zu: %s=%zu, bruteforce=%zu, missing=%zu, extra=%zu\n", + qi, name, got.size(), expected.size(), missing, extra); + } + return all_match; +} + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + float *database = loadRandomData(n_database, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + printf("Loaded %llu database points and %llu query points with dimension %llu\n", + n_database, n_query, dim); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + // Hercules keeps its leaf payloads on disk, so it needs a writable index directory. + daisy::HerculesConfig hercules_config; + hercules_config.index_dir = "/tmp/hercules_range_demo"; + + daisy::Hercules hercules_search(daisy::DistanceType::L2_SQUARED, hercules_config); + hercules_search.setNumThreads(4); + hercules_search.buildIndex(database, n_database, dim); + + std::vector> I; + std::vector> D; + hercules_search.searchIndex(query, n_query, config, I, D); + reportRange("Hercules", r, I, D); + + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.buildIndex(database, n_database, dim); + std::vector> gt_I; + std::vector> gt_D; + ground_truth.searchIndex(query, n_query, config, gt_I, gt_D); + + printf("\nAll queries match: %s\n", + matchesBruteforce("Hercules", I, gt_I) ? "yes" : "no"); + + delete[] database; + delete[] query; + + return 0; +} diff --git a/demos/demo_LbBruteforce_Range.cpp b/demos/demo_LbBruteforce_Range.cpp new file mode 100644 index 0000000..a2fd6d1 --- /dev/null +++ b/demos/demo_LbBruteforce_Range.cpp @@ -0,0 +1,98 @@ +// LbBruteforce range search: report every series within squared-L2 distance r of a query. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +// Range search returns an unordered set of hits, so compare sets rather than positions. +// "missing" are hits brute force found and the index did not; "extra" are the other way +// round, which for an exact index means false positives. +static bool matchesBruteforce(const char *name, + const std::vector> &I, + const std::vector> >_I) +{ + bool all_match = true; + printf("\n=== Verification against brute force ===\n"); + for (size_t qi = 0; qi < I.size(); qi++) + { + const std::set got(I[qi].begin(), I[qi].end()); + const std::set expected(gt_I[qi].begin(), gt_I[qi].end()); + + size_t missing = 0; + for (daisy::idx_t idx : expected) + if (got.count(idx) == 0) missing++; + size_t extra = 0; + for (daisy::idx_t idx : got) + if (expected.count(idx) == 0) extra++; + + all_match = all_match && (missing == 0 && extra == 0); + printf("Query %zu: %s=%zu, bruteforce=%zu, missing=%zu, extra=%zu\n", + qi, name, got.size(), expected.size(), missing, extra); + } + return all_match; +} + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + float *database = loadRandomData(n_database, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + printf("Loaded %llu database points and %llu query points with dimension %llu\n", + n_database, n_query, dim); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + daisy::LbBruteforce lb_search(daisy::DistanceType::L2_SQUARED); + lb_search.setNumThreads(4); + lb_search.buildIndex(database, n_database, dim); + + std::vector> I; + std::vector> D; + lb_search.searchIndex(query, n_query, config, I, D); + reportRange("LbBruteforce", r, I, D); + + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.buildIndex(database, n_database, dim); + std::vector> gt_I; + std::vector> gt_D; + ground_truth.searchIndex(query, n_query, config, gt_I, gt_D); + + printf("\nAll queries match: %s\n", + matchesBruteforce("LbBruteforce", I, gt_I) ? "yes" : "no"); + + delete[] database; + delete[] query; + + return 0; +} diff --git a/demos/demo_Messi_Range.cpp b/demos/demo_Messi_Range.cpp new file mode 100644 index 0000000..7d3e177 --- /dev/null +++ b/demos/demo_Messi_Range.cpp @@ -0,0 +1,98 @@ +// Messi range search: report every series within squared-L2 distance r of a query. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +// Range search returns an unordered set of hits, so compare sets rather than positions. +// "missing" are hits brute force found and the index did not; "extra" are the other way +// round, which for an exact index means false positives. +static bool matchesBruteforce(const char *name, + const std::vector> &I, + const std::vector> >_I) +{ + bool all_match = true; + printf("\n=== Verification against brute force ===\n"); + for (size_t qi = 0; qi < I.size(); qi++) + { + const std::set got(I[qi].begin(), I[qi].end()); + const std::set expected(gt_I[qi].begin(), gt_I[qi].end()); + + size_t missing = 0; + for (daisy::idx_t idx : expected) + if (got.count(idx) == 0) missing++; + size_t extra = 0; + for (daisy::idx_t idx : got) + if (expected.count(idx) == 0) extra++; + + all_match = all_match && (missing == 0 && extra == 0); + printf("Query %zu: %s=%zu, bruteforce=%zu, missing=%zu, extra=%zu\n", + qi, name, got.size(), expected.size(), missing, extra); + } + return all_match; +} + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + float *database = loadRandomData(n_database, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + printf("Loaded %llu database points and %llu query points with dimension %llu\n", + n_database, n_query, dim); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + daisy::Messi messi_search(daisy::DistanceType::L2_SQUARED); + messi_search.setNumThreads(4); + messi_search.buildIndex(database, n_database, dim); + + std::vector> I; + std::vector> D; + messi_search.searchIndex(query, n_query, config, I, D); + reportRange("Messi", r, I, D); + + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.buildIndex(database, n_database, dim); + std::vector> gt_I; + std::vector> gt_D; + ground_truth.searchIndex(query, n_query, config, gt_I, gt_D); + + printf("\nAll queries match: %s\n", + matchesBruteforce("Messi", I, gt_I) ? "yes" : "no"); + + delete[] database; + delete[] query; + + return 0; +} diff --git a/demos/demo_Odyssey_Range.cpp b/demos/demo_Odyssey_Range.cpp new file mode 100644 index 0000000..deb583f --- /dev/null +++ b/demos/demo_Odyssey_Range.cpp @@ -0,0 +1,182 @@ +// Odyssey range search: report every series within squared-L2 distance r of a query. +// Odyssey is MPI-distributed: rank 0 stages the data on disk, every rank indexes and +// searches it, and rank 0 reports the merged result. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#if ODYSSEY_MPI +#include +#endif + +#if defined(__unix__) || defined(__APPLE__) +#include +#include +#endif + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +// Range search returns an unordered set of hits, so compare sets rather than positions. +// "missing" are hits brute force found and the index did not; "extra" are the other way +// round, which for an exact index means false positives. +static bool matchesBruteforce(const char *name, + const std::vector> &I, + const std::vector> >_I) +{ + bool all_match = true; + printf("\n=== Verification against brute force ===\n"); + for (size_t qi = 0; qi < I.size(); qi++) + { + const std::set got(I[qi].begin(), I[qi].end()); + const std::set expected(gt_I[qi].begin(), gt_I[qi].end()); + + size_t missing = 0; + for (daisy::idx_t idx : expected) + if (got.count(idx) == 0) missing++; + size_t extra = 0; + for (daisy::idx_t idx : got) + if (expected.count(idx) == 0) extra++; + + all_match = all_match && (missing == 0 && extra == 0); + printf("Query %zu: %s=%zu, bruteforce=%zu, missing=%zu, extra=%zu\n", + qi, name, got.size(), expected.size(), missing, extra); + } + return all_match; +} + +int main(int argc, char *argv[]) +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + std::string temp_db_file = "odyssey_range_db.bin"; + + daisy::OdysseyConfig odyssey_config; + odyssey_config.search_workers = 2; + odyssey_config.index_threads = 4; + odyssey_config.query_threads = 2; + odyssey_config.leaf_size = 1000; + odyssey_config.paa_segments = 16; + odyssey_config.replication_groups = 0; + + daisy::Odyssey odyssey(odyssey_config, daisy::DistanceType::L2_SQUARED, argc, argv); + int rank = odyssey.getMyRank(); + + float *database = nullptr; + if (rank == 0) + { + remove(temp_db_file.c_str()); + database = loadRandomData(n_database, dim, 100, true); + printf("Loaded %llu database points with dimension %llu\n", n_database, dim); + + FILE *fp = fopen(temp_db_file.c_str(), "wb"); + if (fp == nullptr) + { + fprintf(stderr, "Error: Could not create temporary database file\n"); + delete[] database; + return 1; + } + size_t to_write = static_cast(n_database) * static_cast(dim); + size_t written = fwrite(database, sizeof(float), to_write, fp); + fclose(fp); + if (written != to_write) + { + fprintf(stderr, "Error: wrote only %zu floats (expected %zu)\n", written, to_write); + delete[] database; + return 1; + } + } + +#if ODYSSEY_MPI + MPI_Barrier(MPI_COMM_WORLD); +#endif + + static const int PATH_MAX_MPI = 1024; + char path_buf[PATH_MAX_MPI]; + std::memset(path_buf, 0, PATH_MAX_MPI); + if (rank == 0) + { +#if (defined(__unix__) || defined(__APPLE__)) && defined(PATH_MAX) + char resolved[PATH_MAX]; + if (realpath(temp_db_file.c_str(), resolved) != nullptr) + std::strncpy(path_buf, resolved, PATH_MAX_MPI - 1); + else +#endif + std::strncpy(path_buf, temp_db_file.c_str(), PATH_MAX_MPI - 1); + path_buf[PATH_MAX_MPI - 1] = '\0'; + } +#if ODYSSEY_MPI + MPI_Bcast(path_buf, PATH_MAX_MPI, MPI_CHAR, 0, MPI_COMM_WORLD); +#endif + std::string path_to_use(path_buf); + + float *query = loadRandomData(n_query, dim, 50, true); + if (rank == 0) + printf("Loaded %llu query points with dimension %llu\n", n_query, dim); + + daisy::FileDataSource data_source(path_to_use.c_str(), dim, n_database); + odyssey.buildIndex(&data_source); + if (rank == 0) + printf(">>> Finished indexing\n"); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + std::vector> I; + std::vector> D; + odyssey.searchIndex(query, n_query, config, I, D); + + if (rank == 0) + { + reportRange("Odyssey", r, I, D); + + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.buildIndex(database, n_database, dim); + std::vector> gt_I; + std::vector> gt_D; + ground_truth.searchIndex(query, n_query, config, gt_I, gt_D); + + printf("\nAll queries match: %s\n", + matchesBruteforce("Odyssey", I, gt_I) ? "yes" : "no"); + } + + delete[] database; + delete[] query; + if (rank == 0) + remove(path_to_use.c_str()); + +#if ODYSSEY_MPI + MPI_Finalize(); +#endif + + return 0; +} diff --git a/demos/demo_ParIS_Range.cpp b/demos/demo_ParIS_Range.cpp new file mode 100644 index 0000000..a5a011f --- /dev/null +++ b/demos/demo_ParIS_Range.cpp @@ -0,0 +1,112 @@ +// ParIS range search: report every series within squared-L2 distance r of a query. +// ParIS indexes from a file, so the generated data is staged on disk first. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include +#include + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +// Range search returns an unordered set of hits, so compare sets rather than positions. +// "missing" are hits brute force found and the index did not; "extra" are the other way +// round, which for an exact index means false positives. +static bool matchesBruteforce(const char *name, + const std::vector> &I, + const std::vector> >_I) +{ + bool all_match = true; + printf("\n=== Verification against brute force ===\n"); + for (size_t qi = 0; qi < I.size(); qi++) + { + const std::set got(I[qi].begin(), I[qi].end()); + const std::set expected(gt_I[qi].begin(), gt_I[qi].end()); + + size_t missing = 0; + for (daisy::idx_t idx : expected) + if (got.count(idx) == 0) missing++; + size_t extra = 0; + for (daisy::idx_t idx : got) + if (expected.count(idx) == 0) extra++; + + all_match = all_match && (missing == 0 && extra == 0); + printf("Query %zu: %s=%zu, bruteforce=%zu, missing=%zu, extra=%zu\n", + qi, name, got.size(), expected.size(), missing, extra); + } + return all_match; +} + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + float *database = loadRandomData(n_database, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + printf("Loaded %llu database points and %llu query points with dimension %llu\n", + n_database, n_query, dim); + + std::string temp_db_file = "/tmp/paris_range_db.bin"; + FILE *fp = fopen(temp_db_file.c_str(), "wb"); + if (fp == nullptr) { + fprintf(stderr, "Error: Could not create temporary database file\n"); + delete[] database; + delete[] query; + return 1; + } + fwrite(database, sizeof(float), n_database * dim, fp); + fclose(fp); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + daisy::ParIS paris_search(daisy::DistanceType::L2_SQUARED); + paris_search.setNumThreads(4); + paris_search.buildIndex(temp_db_file, dim, n_database); + + std::vector> I; + std::vector> D; + paris_search.searchIndex(query, n_query, config, I, D); + reportRange("ParIS", r, I, D); + + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.buildIndex(database, n_database, dim); + std::vector> gt_I; + std::vector> gt_D; + ground_truth.searchIndex(query, n_query, config, gt_I, gt_D); + + printf("\nAll queries match: %s\n", + matchesBruteforce("ParIS", I, gt_I) ? "yes" : "no"); + + delete[] database; + delete[] query; + remove(temp_db_file.c_str()); + + return 0; +} diff --git a/demos/demo_Sing_Range.cpp b/demos/demo_Sing_Range.cpp new file mode 100644 index 0000000..2f7c067 --- /dev/null +++ b/demos/demo_Sing_Range.cpp @@ -0,0 +1,110 @@ +// Sing range search: report every series within squared-L2 distance r of a query. +// Sing is the CUDA-accelerated index, so it is only built when CUDA is available. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include +#include + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +// Range search returns an unordered set of hits, so compare sets rather than positions. +// "missing" are hits brute force found and the index did not; "extra" are the other way +// round, which for an exact index means false positives. +static bool matchesBruteforce(const char *name, + const std::vector> &I, + const std::vector> >_I) +{ + bool all_match = true; + printf("\n=== Verification against brute force ===\n"); + for (size_t qi = 0; qi < I.size(); qi++) + { + const std::set got(I[qi].begin(), I[qi].end()); + const std::set expected(gt_I[qi].begin(), gt_I[qi].end()); + + size_t missing = 0; + for (daisy::idx_t idx : expected) + if (got.count(idx) == 0) missing++; + size_t extra = 0; + for (daisy::idx_t idx : got) + if (expected.count(idx) == 0) extra++; + + all_match = all_match && (missing == 0 && extra == 0); + printf("Query %zu: %s=%zu, bruteforce=%zu, missing=%zu, extra=%zu\n", + qi, name, got.size(), expected.size(), missing, extra); + } + return all_match; +} + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + float *database = loadRandomData(n_database, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + printf("Loaded %llu database points and %llu query points with dimension %llu\n", + n_database, n_query, dim); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + daisy::Sing sing_search(daisy::DistanceType::L2_SQUARED); + + daisy::InMemoryDataSource data_source(database, n_database, dim); + auto t0 = std::chrono::steady_clock::now(); + sing_search.buildIndex(&data_source); + auto t1 = std::chrono::steady_clock::now(); + double build_ms = 1e-6 * (double)std::chrono::duration_cast(t1 - t0).count(); + printf("buildIndex done in %.2f ms\n", build_ms); + + std::vector> I; + std::vector> D; + auto t_search0 = std::chrono::steady_clock::now(); + sing_search.searchIndex(query, n_query, config, I, D); + auto t_search1 = std::chrono::steady_clock::now(); + double search_ms = 1e-6 * (double)std::chrono::duration_cast(t_search1 - t_search0).count(); + printf("Range search done in %.2f ms (%.2f ms/query)\n", search_ms, search_ms / (double)n_query); + + reportRange("Sing", r, I, D); + + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.buildIndex(database, n_database, dim); + std::vector> gt_I; + std::vector> gt_D; + ground_truth.searchIndex(query, n_query, config, gt_I, gt_D); + + printf("\nAll queries match: %s\n", + matchesBruteforce("Sing", I, gt_I) ? "yes" : "no"); + + delete[] database; + delete[] query; + + return 0; +} diff --git a/demos/demo_Sofa_Range.cpp b/demos/demo_Sofa_Range.cpp new file mode 100644 index 0000000..7a76f18 --- /dev/null +++ b/demos/demo_Sofa_Range.cpp @@ -0,0 +1,98 @@ +// Sofa range search: report every series within squared-L2 distance r of a query. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +#include +#include +#include +#include + +static void reportRange(const char *name, float r, + const std::vector> &I, + const std::vector> &D) +{ + printf("=== %s range search (squared L2 <= %.1f) ===\n", name, r); + for (size_t qi = 0; qi < I.size(); qi++) + { + printf("Query %zu: %zu hits", qi, I[qi].size()); + if (!I[qi].empty()) + { + const auto minmax = std::minmax_element(D[qi].begin(), D[qi].end()); + printf(" [closest=%.4f, farthest=%.4f]", *minmax.first, *minmax.second); + } + printf("\n"); + } +} + +// Range search returns an unordered set of hits, so compare sets rather than positions. +// "missing" are hits brute force found and the index did not; "extra" are the other way +// round, which for an exact index means false positives. +static bool matchesBruteforce(const char *name, + const std::vector> &I, + const std::vector> >_I) +{ + bool all_match = true; + printf("\n=== Verification against brute force ===\n"); + for (size_t qi = 0; qi < I.size(); qi++) + { + const std::set got(I[qi].begin(), I[qi].end()); + const std::set expected(gt_I[qi].begin(), gt_I[qi].end()); + + size_t missing = 0; + for (daisy::idx_t idx : expected) + if (got.count(idx) == 0) missing++; + size_t extra = 0; + for (daisy::idx_t idx : got) + if (expected.count(idx) == 0) extra++; + + all_match = all_match && (missing == 0 && extra == 0); + printf("Query %zu: %s=%zu, bruteforce=%zu, missing=%zu, extra=%zu\n", + qi, name, got.size(), expected.size(), missing, extra); + } + return all_match; +} + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + + // Independent z-normalized series sit around 2*dim apart, so a radius below that + // keeps the result sets small without leaving every query empty. + const float r = 0.60f * 2.0f * dim; + + float *database = loadRandomData(n_database, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + printf("Loaded %llu database points and %llu query points with dimension %llu\n", + n_database, n_query, dim); + + daisy::SearchConfig config; + config.type = daisy::QueryType::RANGE; + config.r = r; + + daisy::Sofa sofa_search(daisy::DistanceType::L2_SQUARED); + sofa_search.setNumThreads(4); + sofa_search.buildIndex(database, n_database, dim); + + std::vector> I; + std::vector> D; + sofa_search.searchIndex(query, n_query, config, I, D); + reportRange("Sofa", r, I, D); + + daisy::BruteForceSearch ground_truth(daisy::DistanceType::L2_SQUARED); + ground_truth.buildIndex(database, n_database, dim); + std::vector> gt_I; + std::vector> gt_D; + ground_truth.searchIndex(query, n_query, config, gt_I, gt_D); + + printf("\nAll queries match: %s\n", + matchesBruteforce("Sofa", I, gt_I) ? "yes" : "no"); + + delete[] database; + delete[] query; + + return 0; +} diff --git a/docs/demos-guide.md b/docs/demos-guide.md index 6b80be5..ca45601 100644 --- a/docs/demos-guide.md +++ b/docs/demos-guide.md @@ -9,6 +9,14 @@ Most demos follow the same batch pattern: `buildIndex(...)` once, then `searchIn `demo_LbBruteforce_Streaming`, `demo_Messi_Streaming`, and `demo_Coconut_Streaming` for live-index examples. +Range search is exposed through `SearchConfig` and is implemented by **Bruteforce**, +**LbBruteforce**, **MESSI**, **Coconut**, **ParIS**, **Fresh**, **DumpyOS**, **Hercules**, +**Sofa**, **Sing**, and **Odyssey**. Each one has a `demo__Range` C++ demo: set +`config.type = daisy::QueryType::RANGE` and `config.r`, then read the per-query hit lists from +the `std::vector>` overload of `searchIndex(...)`. Unlike top-k, a range query +returns a variable number of unordered hits per query, so the demos compare result *sets* +against brute force rather than positions. + ## Demo Program Structure All demos are located in the [`demos/`](../demos/) directory. @@ -21,7 +29,7 @@ Each demo can be customized by modifying: - **Dataset Size**: `n_database`, `n_query` - **Dimensionality**: `dim` (time series length) -- **Search Parameters**: `k` (number of neighbors) +- **Search Parameters**: `k` (number of neighbors), `r` (range-query radius) - **Algorithm Parameters**: Thread count, distance metrics, etc. ### Data Sources From ef797b96e73f15932b0aafb3c08a056b1ee397cd Mon Sep 17 00:00:00 2001 From: sophisid Date: Thu, 17 Sep 2026 13:16:12 +0200 Subject: [PATCH 4/4] fix breakpoints issue --- lib/algos/LbBruteforce.cpp | 12 +- lib/algos/Messi.cpp | 23 +-- lib/distance_computers/DistanceComputer.cpp | 36 ++-- lib/distance_computers/DistanceComputer.hpp | 18 +- lib/isax/SAX.cpp | 196 +++++++++++--------- lib/isax/SAX.hpp | 24 ++- tests/test_equidepth.cpp | 118 ++++++++++++ 7 files changed, 298 insertions(+), 129 deletions(-) diff --git a/lib/algos/LbBruteforce.cpp b/lib/algos/LbBruteforce.cpp index 80359f1..ac211ee 100644 --- a/lib/algos/LbBruteforce.cpp +++ b/lib/algos/LbBruteforce.cpp @@ -272,7 +272,6 @@ namespace daisy { if (!validateSearchParams(k, n_query)) return; - activateBreakpoints(); #pragma omp parallel num_threads(num_threads) { @@ -308,7 +307,8 @@ namespace daisy index->settings->paa_segments, MINVAL, MAXVAL, - index->settings->mindist_sqrt); + index->settings->mindist_sqrt, + index->settings->breakpoints); if (minimum_distance < bound) { float dist = this->distance_computer->compute_dist_SIMD(const_cast(q_vec), @@ -345,7 +345,6 @@ namespace daisy { if (!validateSearchParams(k, n_query)) return; - activateBreakpoints(); #pragma omp parallel num_threads(num_threads) { @@ -401,7 +400,8 @@ namespace daisy index->settings->paa_segments, MINVAL, MAXVAL, - index->settings->mindist_sqrt); + index->settings->mindist_sqrt, + index->settings->breakpoints); if (minimum_distance < bound) { @@ -452,7 +452,6 @@ namespace daisy throw std::runtime_error("LbBruteforce index must be built before searching"); if (n_query == 0) throw std::invalid_argument("n_query must be greater than 0"); - activateBreakpoints(); float r = config.r; I.assign(n_query, {}); @@ -486,7 +485,8 @@ namespace daisy index->settings->paa_segments, MINVAL, MAXVAL, - index->settings->mindist_sqrt); + index->settings->mindist_sqrt, + index->settings->breakpoints); if (minimum_distance <= r) { float dist = this->distance_computer->compute_dist_SIMD( diff --git a/lib/algos/Messi.cpp b/lib/algos/Messi.cpp index 4134c71..14e7cd3 100644 --- a/lib/algos/Messi.cpp +++ b/lib/algos/Messi.cpp @@ -33,7 +33,8 @@ namespace daisy index->settings->sax_alphabet_cardinality, index->settings->paa_segments, MINVAL, MAXVAL, - index->settings->mindist_sqrt); + index->settings->mindist_sqrt, + index->settings->breakpoints); if (distance < bsf) { @@ -80,7 +81,8 @@ namespace daisy index->settings->sax_alphabet_cardinality, index->settings->paa_segments, MINVAL, MAXVAL, - index->settings->mindist_sqrt); + index->settings->mindist_sqrt, + index->settings->breakpoints); if (distance <= bsf) { if (node->is_leaf) @@ -165,7 +167,8 @@ namespace daisy index->settings->sax_bit_cardinality, index->settings->sax_alphabet_cardinality, index->settings->paa_segments, MINVAL, MAXVAL, - index->settings->mindist_sqrt); + index->settings->mindist_sqrt, + index->settings->breakpoints); if (distmin <= pq_bsf->knn[pq_bsf->k - 1]) { float dist = ts_euclidean_distance_SIMD(query, &(rawfile[*node->buffer->partial_position_buffer[i]]), @@ -210,7 +213,8 @@ namespace daisy index->settings->sax_bit_cardinality, index->settings->sax_alphabet_cardinality, index->settings->paa_segments, MINVAL, MAXVAL, - index->settings->mindist_sqrt); + index->settings->mindist_sqrt, + index->settings->breakpoints); if (distmin <= bsf) { @@ -531,7 +535,8 @@ namespace daisy index->settings->sax_bit_cardinality, index->settings->sax_alphabet_cardinality, index->settings->paa_segments, MINVAL, MAXVAL, - index->settings->mindist_sqrt); + index->settings->mindist_sqrt, + index->settings->breakpoints); if (distmin <= r) { float dist = ts_euclidean_distance_SIMD(query, &(rawfile[*node->buffer->partial_position_buffer[i]]), @@ -1088,11 +1093,6 @@ namespace daisy void Messi::searchIndexL2Squared(const float *query, const idx_t n_query, const idx_t k, idx_t *I, float *D) { - // Lower bounds read the active breakpoints, which are global: interleaving - // searches with inserts (or with another live index) can leave a different - // table installed, so reinstall ours before every search. - activateBreakpoints(); - ts_type *paa = (ts_type *)malloc(sizeof(ts_type) * index->settings->paa_segments); node_list nodelist; @@ -1167,8 +1167,6 @@ namespace daisy void Messi::searchIndexDTW(const float *query, const idx_t n_query, const idx_t k, idx_t *I, float *D) { - activateBreakpoints(); - isax_index *index = this->index; node_list nodelist; @@ -1308,7 +1306,6 @@ namespace daisy return; } - activateBreakpoints(); ts_type *paa = (ts_type *)malloc(sizeof(ts_type) * index->settings->paa_segments); node_list nodelist; nodelist.nlist = (isax_node **)malloc(sizeof(isax_node *) * (int)pow(2, index->settings->paa_segments)); diff --git a/lib/distance_computers/DistanceComputer.cpp b/lib/distance_computers/DistanceComputer.cpp index 0ed38df..def2c43 100644 --- a/lib/distance_computers/DistanceComputer.cpp +++ b/lib/distance_computers/DistanceComputer.cpp @@ -89,7 +89,8 @@ namespace daisy int paa_segments, float minval, float maxval, - bool mindist_sqrt) + bool mindist_sqrt, + const float *bp) { return minidist_paa_to_isax_rawa_SIMD( const_cast(q_paa), @@ -100,7 +101,8 @@ namespace daisy paa_segments, minval, maxval, - mindist_sqrt); + mindist_sqrt, + bp); } void DistanceComputer::compute_paa_from_ts(const float *ts, @@ -149,12 +151,14 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp) { return minidist_paa_to_isax(paa, sax, sax_cardinalities, max_bit_cardinality, max_cardinality, number_of_segments, min_val, max_val, - ratio_sqrt); + ratio_sqrt, + bp); } float DistanceComputer::wrap_minidist_paa_to_isax_raw_SIMD(float *paa, sax_type *sax, @@ -164,12 +168,14 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp) { return minidist_paa_to_isax_raw_SIMD(paa, sax, sax_cardinalities, max_bit_cardinality, max_cardinality, number_of_segments, min_val, max_val, - ratio_sqrt); + ratio_sqrt, + bp); } float DistanceComputer::wrap_ts_euclidean_distance(ts_type *t, ts_type *s, int size, float bound) @@ -189,12 +195,14 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp) { return minidist_paa_to_isax_rawa_SIMD(paa, sax, sax_cardinalities, max_bit_cardinality, max_cardinality, number_of_segments, min_val, max_val, - ratio_sqrt); + ratio_sqrt, + bp); } float DistanceComputer::wrap_minidist_paa_to_isax_raw_DTW_SIMD(float *paaU, float *paaL, sax_type *sax, @@ -204,12 +212,14 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp) { return minidist_paa_to_isax_raw_DTW_SIMD(paaU, paaL, sax, sax_cardinalities, max_bit_cardinality, max_cardinality, number_of_segments, min_val, max_val, - ratio_sqrt); + ratio_sqrt, + bp); } float DistanceComputer::wrap_lb_keogh_data_bound(float *qo, float *tu, float *tl, float *cb, int len, float bsf) @@ -229,12 +239,14 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp) { return minidist_paa_to_isax_DTW(paaU, paaL, sax, sax_cardinalities, max_bit_cardinality, max_cardinality, number_of_segments, min_val, max_val, - ratio_sqrt); + ratio_sqrt, + bp); } // DTW distance methods implementation diff --git a/lib/distance_computers/DistanceComputer.hpp b/lib/distance_computers/DistanceComputer.hpp index f5ab918..398a70a 100644 --- a/lib/distance_computers/DistanceComputer.hpp +++ b/lib/distance_computers/DistanceComputer.hpp @@ -90,7 +90,8 @@ namespace daisy int paa_segments, float minval, float maxval, - bool mindist_sqrt); + bool mindist_sqrt, + const float *bp = nullptr); void compute_paa_from_ts(const float *ts, ts_type *paa, @@ -117,7 +118,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp = nullptr); float wrap_minidist_paa_to_isax_raw_SIMD(float *paa, sax_type *sax, sax_type *sax_cardinalities, @@ -126,7 +128,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp = nullptr); float wrap_ts_euclidean_distance(ts_type *t, ts_type *s, int size, float bound); @@ -139,7 +142,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp = nullptr); float wrap_minidist_paa_to_isax_raw_DTW_SIMD(float *paaU, float *paaL, sax_type *sax, sax_type *sax_cardinalities, @@ -148,7 +152,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp = nullptr); float wrap_lb_keogh_data_bound(float *qo, float *tu, float *tl, float *cb, int len, float bsf); @@ -161,7 +166,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp = nullptr); }; } // namespace daisy diff --git a/lib/isax/SAX.cpp b/lib/isax/SAX.cpp index 2116c06..bc537bc 100644 --- a/lib/isax/SAX.cpp +++ b/lib/isax/SAX.cpp @@ -153,8 +153,12 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp) { + // A null table means the caller has not been migrated off the global yet. + if (bp == nullptr) + bp = daisy_active_breakpoints; float distance = 0; // TODO: Store offset in index settings. and pass index settings as parameter. @@ -182,7 +186,7 @@ namespace daisy } else { - breakpoint_lower = daisy_active_breakpoints[offset + region_lower - 1]; + breakpoint_lower = bp[offset + region_lower - 1]; } if (region_upper == max_cardinality - 1) { @@ -190,7 +194,7 @@ namespace daisy } else { - breakpoint_upper = daisy_active_breakpoints[offset + region_upper]; + breakpoint_upper = bp[offset + region_upper]; } if (breakpoint_lower > paa[i]) @@ -215,8 +219,12 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp) { + // A null table means the caller has not been migrated off the global yet. + if (bp == nullptr) + bp = daisy_active_breakpoints; int region_upper[16], region_lower[16]; float distancef[16]; @@ -277,43 +285,43 @@ namespace daisy __m256 minvalv = _mm256_set1_ps(min_val); //__m256 lsax_breakpoints_shiftv_0 _mm256_i32gather_ps (sax_breakpoints, __m256i vindex, const int scale) - __m256 lsax_breakpoints_shiftv_0 = _mm256_set_ps(daisy_active_breakpoints[offset + region_lower[7] - 1], - daisy_active_breakpoints[offset + region_lower[6] - 1], - daisy_active_breakpoints[offset + region_lower[5] - 1], - daisy_active_breakpoints[offset + region_lower[4] - 1], - daisy_active_breakpoints[offset + region_lower[3] - 1], - daisy_active_breakpoints[offset + region_lower[2] - 1], - daisy_active_breakpoints[offset + region_lower[1] - 1], - daisy_active_breakpoints[offset + region_lower[0] - 1]); - __m256 lsax_breakpoints_shiftv_1 = _mm256_set_ps(daisy_active_breakpoints[offset + region_lower[15] - 1], - daisy_active_breakpoints[offset + region_lower[14] - 1], - daisy_active_breakpoints[offset + region_lower[13] - 1], - daisy_active_breakpoints[offset + region_lower[12] - 1], - daisy_active_breakpoints[offset + region_lower[11] - 1], - daisy_active_breakpoints[offset + region_lower[10] - 1], - daisy_active_breakpoints[offset + region_lower[9] - 1], - daisy_active_breakpoints[offset + region_lower[8] - 1]); + __m256 lsax_breakpoints_shiftv_0 = _mm256_set_ps(bp[offset + region_lower[7] - 1], + bp[offset + region_lower[6] - 1], + bp[offset + region_lower[5] - 1], + bp[offset + region_lower[4] - 1], + bp[offset + region_lower[3] - 1], + bp[offset + region_lower[2] - 1], + bp[offset + region_lower[1] - 1], + bp[offset + region_lower[0] - 1]); + __m256 lsax_breakpoints_shiftv_1 = _mm256_set_ps(bp[offset + region_lower[15] - 1], + bp[offset + region_lower[14] - 1], + bp[offset + region_lower[13] - 1], + bp[offset + region_lower[12] - 1], + bp[offset + region_lower[11] - 1], + bp[offset + region_lower[10] - 1], + bp[offset + region_lower[9] - 1], + bp[offset + region_lower[8] - 1]); __m256 breakpoint_lowerv_0 = (__m256)_mm256_or_si256(_mm256_and_si256(lower_juge_zerov_0, (__m256i)minvalv), _mm256_and_si256(lower_juge_nzerov_0, (__m256i)lsax_breakpoints_shiftv_0)); __m256 breakpoint_lowerv_1 = (__m256)_mm256_or_si256(_mm256_and_si256(lower_juge_zerov_1, (__m256i)minvalv), _mm256_and_si256(lower_juge_nzerov_1, (__m256i)lsax_breakpoints_shiftv_1)); // uper - __m256 usax_breakpoints_shiftv_0 = _mm256_set_ps(daisy_active_breakpoints[offset + region_upper[7]], - daisy_active_breakpoints[offset + region_upper[6]], - daisy_active_breakpoints[offset + region_upper[5]], - daisy_active_breakpoints[offset + region_upper[4]], - daisy_active_breakpoints[offset + region_upper[3]], - daisy_active_breakpoints[offset + region_upper[2]], - daisy_active_breakpoints[offset + region_upper[1]], - daisy_active_breakpoints[offset + region_upper[0]]); - __m256 usax_breakpoints_shiftv_1 = _mm256_set_ps(daisy_active_breakpoints[offset + region_upper[15]], - daisy_active_breakpoints[offset + region_upper[14]], - daisy_active_breakpoints[offset + region_upper[13]], - daisy_active_breakpoints[offset + region_upper[12]], - daisy_active_breakpoints[offset + region_upper[11]], - daisy_active_breakpoints[offset + region_upper[10]], - daisy_active_breakpoints[offset + region_upper[9]], - daisy_active_breakpoints[offset + region_upper[8]]); + __m256 usax_breakpoints_shiftv_0 = _mm256_set_ps(bp[offset + region_upper[7]], + bp[offset + region_upper[6]], + bp[offset + region_upper[5]], + bp[offset + region_upper[4]], + bp[offset + region_upper[3]], + bp[offset + region_upper[2]], + bp[offset + region_upper[1]], + bp[offset + region_upper[0]]); + __m256 usax_breakpoints_shiftv_1 = _mm256_set_ps(bp[offset + region_upper[15]], + bp[offset + region_upper[14]], + bp[offset + region_upper[13]], + bp[offset + region_upper[12]], + bp[offset + region_upper[11]], + bp[offset + region_upper[10]], + bp[offset + region_upper[9]], + bp[offset + region_upper[8]]); __m256i upper_juge_maxv_0 = _mm256_cmpeq_epi32(region_upperv_0, _mm256_set1_epi32(max_cardinality - 1)); __m256i upper_juge_maxv_1 = _mm256_cmpeq_epi32(region_upperv_1, _mm256_set1_epi32(max_cardinality - 1)); @@ -378,8 +386,12 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp_max) { + // A null table means the caller has not been migrated off the global yet. + if (bp_max == nullptr) + bp_max = daisy_active_breakpoints_max; int region_upper[16], region_lower[16]; float distancef[16]; int offset = 0; @@ -439,8 +451,8 @@ namespace daisy __m256 minvalv = _mm256_set1_ps((float)min_val); - __m256 lsax_breakpoints_shiftv_0 = _mm256_i32gather_ps(daisy_active_breakpoints_max, region_lowerv_0, 4); - __m256 lsax_breakpoints_shiftv_1 = _mm256_i32gather_ps(daisy_active_breakpoints_max, region_lowerv_1, 4); + __m256 lsax_breakpoints_shiftv_0 = _mm256_i32gather_ps(bp_max, region_lowerv_0, 4); + __m256 lsax_breakpoints_shiftv_1 = _mm256_i32gather_ps(bp_max, region_lowerv_1, 4); __m256 breakpoint_lowerv_0 = (__m256)_mm256_or_si256( _mm256_and_si256(lower_juge_zerov_0, (__m256i)minvalv), @@ -449,8 +461,8 @@ namespace daisy _mm256_and_si256(lower_juge_zerov_1, (__m256i)minvalv), _mm256_and_si256(lower_juge_nzerov_1, (__m256i)lsax_breakpoints_shiftv_1)); - __m256 usax_breakpoints_shiftv_0 = _mm256_i32gather_ps(daisy_active_breakpoints_max, region_upperv_0, 4); - __m256 usax_breakpoints_shiftv_1 = _mm256_i32gather_ps(daisy_active_breakpoints_max, region_upperv_1, 4); + __m256 usax_breakpoints_shiftv_0 = _mm256_i32gather_ps(bp_max, region_upperv_0, 4); + __m256 usax_breakpoints_shiftv_1 = _mm256_i32gather_ps(bp_max, region_upperv_1, 4); __m256i upper_juge_maxv_0 = _mm256_cmpeq_epi32(region_upperv_0, _mm256_set1_epi32(max_cardinality - 1)); __m256i upper_juge_maxv_1 = _mm256_cmpeq_epi32(region_upperv_1, _mm256_set1_epi32(max_cardinality - 1)); @@ -582,8 +594,12 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp) { + // A null table means the caller has not been migrated off the global yet. + if (bp == nullptr) + bp = daisy_active_breakpoints; int region_upper[16], region_lower[16]; float distancef[16]; @@ -660,15 +676,15 @@ namespace daisy __m256 minvalv = _mm256_set1_ps(min_val); - __m256 lsax_breakpoints_shiftv_0 = _mm256_i32gather_ps(daisy_active_breakpoints, region_lowerv_0_offset, 4); - __m256 lsax_breakpoints_shiftv_1 = _mm256_i32gather_ps(daisy_active_breakpoints, region_lowerv_1_offset, 4); + __m256 lsax_breakpoints_shiftv_0 = _mm256_i32gather_ps(bp, region_lowerv_0_offset, 4); + __m256 lsax_breakpoints_shiftv_1 = _mm256_i32gather_ps(bp, region_lowerv_1_offset, 4); __m256 breakpoint_lowerv_0 = (__m256)_mm256_or_si256(_mm256_and_si256(lower_juge_zerov_0, (__m256i)minvalv), _mm256_and_si256(lower_juge_nzerov_0, (__m256i)lsax_breakpoints_shiftv_0)); __m256 breakpoint_lowerv_1 = (__m256)_mm256_or_si256(_mm256_and_si256(lower_juge_zerov_1, (__m256i)minvalv), _mm256_and_si256(lower_juge_nzerov_1, (__m256i)lsax_breakpoints_shiftv_1)); // uper - __m256 usax_breakpoints_shiftv_0 = _mm256_i32gather_ps(daisy_active_breakpoints, region_upperv_0_offset, 4); - __m256 usax_breakpoints_shiftv_1 = _mm256_i32gather_ps(daisy_active_breakpoints, region_upperv_1_offset, 4); + __m256 usax_breakpoints_shiftv_0 = _mm256_i32gather_ps(bp, region_upperv_0_offset, 4); + __m256 usax_breakpoints_shiftv_1 = _mm256_i32gather_ps(bp, region_upperv_1_offset, 4); __m256i upper_juge_maxv_0 = _mm256_cmpeq_epi32(region_upperv_0, _mm256_set1_epi32(max_cardinality - 1)); __m256i upper_juge_maxv_1 = _mm256_cmpeq_epi32(region_upperv_1, _mm256_set1_epi32(max_cardinality - 1)); @@ -733,8 +749,12 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp) { + // A null table means the caller has not been migrated off the global yet. + if (bp == nullptr) + bp = daisy_active_breakpoints; int region_upper[16], region_lower[16]; float distancef[16]; @@ -804,47 +824,47 @@ namespace daisy __m256 minvalv = _mm256_set1_ps(min_val); - __m256 lsax_breakpoints_shiftv_0 = _mm256_i32gather_ps(daisy_active_breakpoints, region_lowerv_0_offset, 4); - //__m256 lsax_breakpoints_shiftv_0= _mm256_set_ps (daisy_active_breakpoints[region_lower[7]], - // daisy_active_breakpoints[region_lower[6]], - // daisy_active_breakpoints[region_lower[5]], - // daisy_active_breakpoints[region_lower[4]], - // daisy_active_breakpoints[region_lower[3]], - // daisy_active_breakpoints[region_lower[2]], - // daisy_active_breakpoints[region_lower[1]], - // daisy_active_breakpoints[region_lower[0]]); - __m256 lsax_breakpoints_shiftv_1 = _mm256_i32gather_ps(daisy_active_breakpoints, region_lowerv_1_offset, 4); - //__m256 lsax_breakpoints_shiftv_1= _mm256_set_ps (daisy_active_breakpoints[region_lower[15]], - // daisy_active_breakpoints[region_lower[14]], - // daisy_active_breakpoints[region_lower[13]], - // daisy_active_breakpoints[region_lower[12]], - // daisy_active_breakpoints[region_lower[11]], - // daisy_active_breakpoints[region_lower[10]], - // daisy_active_breakpoints[region_lower[9]], - // daisy_active_breakpoints[region_lower[8]]); + __m256 lsax_breakpoints_shiftv_0 = _mm256_i32gather_ps(bp, region_lowerv_0_offset, 4); + //__m256 lsax_breakpoints_shiftv_0= _mm256_set_ps (bp[region_lower[7]], + // bp[region_lower[6]], + // bp[region_lower[5]], + // bp[region_lower[4]], + // bp[region_lower[3]], + // bp[region_lower[2]], + // bp[region_lower[1]], + // bp[region_lower[0]]); + __m256 lsax_breakpoints_shiftv_1 = _mm256_i32gather_ps(bp, region_lowerv_1_offset, 4); + //__m256 lsax_breakpoints_shiftv_1= _mm256_set_ps (bp[region_lower[15]], + // bp[region_lower[14]], + // bp[region_lower[13]], + // bp[region_lower[12]], + // bp[region_lower[11]], + // bp[region_lower[10]], + // bp[region_lower[9]], + // bp[region_lower[8]]); __m256 breakpoint_lowerv_0 = (__m256)_mm256_or_si256(_mm256_and_si256(lower_juge_zerov_0, (__m256i)minvalv), _mm256_and_si256(lower_juge_nzerov_0, (__m256i)lsax_breakpoints_shiftv_0)); __m256 breakpoint_lowerv_1 = (__m256)_mm256_or_si256(_mm256_and_si256(lower_juge_zerov_1, (__m256i)minvalv), _mm256_and_si256(lower_juge_nzerov_1, (__m256i)lsax_breakpoints_shiftv_1)); // uper - __m256 usax_breakpoints_shiftv_0 = _mm256_i32gather_ps(daisy_active_breakpoints, region_upperv_0_offset, 4); - //__m256 usax_breakpoints_shiftv_0= _mm256_set_ps (daisy_active_breakpoints[region_upper[7]], - // daisy_active_breakpoints[region_upper[6]], - // daisy_active_breakpoints[region_upper[5]], - // daisy_active_breakpoints[region_upper[4]], - // daisy_active_breakpoints[region_upper[3]], - // daisy_active_breakpoints[region_upper[2]], - // daisy_active_breakpoints[region_upper[1]], - // daisy_active_breakpoints[region_upper[0]]); - __m256 usax_breakpoints_shiftv_1 = _mm256_i32gather_ps(daisy_active_breakpoints, region_upperv_1_offset, 4); - //__m256 usax_breakpoints_shiftv_1= _mm256_set_ps (daisy_active_breakpoints[region_upper[15]], - // daisy_active_breakpoints[region_upper[14]], - // daisy_active_breakpoints[region_upper[13]], - // daisy_active_breakpoints[region_upper[12]], - // daisy_active_breakpoints[region_upper[11]], - // daisy_active_breakpoints[region_upper[10]], - // daisy_active_breakpoints[region_upper[9]], - // daisy_active_breakpoints[region_upper[8]]); + __m256 usax_breakpoints_shiftv_0 = _mm256_i32gather_ps(bp, region_upperv_0_offset, 4); + //__m256 usax_breakpoints_shiftv_0= _mm256_set_ps (bp[region_upper[7]], + // bp[region_upper[6]], + // bp[region_upper[5]], + // bp[region_upper[4]], + // bp[region_upper[3]], + // bp[region_upper[2]], + // bp[region_upper[1]], + // bp[region_upper[0]]); + __m256 usax_breakpoints_shiftv_1 = _mm256_i32gather_ps(bp, region_upperv_1_offset, 4); + //__m256 usax_breakpoints_shiftv_1= _mm256_set_ps (bp[region_upper[15]], + // bp[region_upper[14]], + // bp[region_upper[13]], + // bp[region_upper[12]], + // bp[region_upper[11]], + // bp[region_upper[10]], + // bp[region_upper[9]], + // bp[region_upper[8]]); __m256i upper_juge_maxv_0 = _mm256_cmpeq_epi32(region_upperv_0, _mm256_set1_epi32(max_cardinality - 1)); __m256i upper_juge_maxv_1 = _mm256_cmpeq_epi32(region_upperv_1, _mm256_set1_epi32(max_cardinality - 1)); @@ -1090,8 +1110,12 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt) + float ratio_sqrt, + const float *bp) { + // A null table means the caller has not been migrated off the global yet. + if (bp == nullptr) + bp = daisy_active_breakpoints; float distance = 0; // TODO: Store offset in index settings. and pass index settings as parameter. @@ -1121,7 +1145,7 @@ namespace daisy } else { - breakpoint_lower = daisy_active_breakpoints[offset + region_lower - 1]; + breakpoint_lower = bp[offset + region_lower - 1]; } if (region_upper == max_cardinality - 1) { @@ -1129,7 +1153,7 @@ namespace daisy } else { - breakpoint_upper = daisy_active_breakpoints[offset + region_upper]; + breakpoint_upper = bp[offset + region_upper]; } if (breakpoint_lower > paaU[i]) diff --git a/lib/isax/SAX.hpp b/lib/isax/SAX.hpp index 7af0984..4cea42b 100644 --- a/lib/isax/SAX.hpp +++ b/lib/isax/SAX.hpp @@ -11,6 +11,12 @@ namespace daisy { // Process-global active breakpoints (triangular / flat-max). Default to the Gaussian // tables; an index installs its own via set_active_breakpoints() before build/search. + // + // Being process-global, these are only correct while a single index is live: two indices + // with different (equi-depth) tables would read each other's. They are being retired. The + // lower-bound helpers below now take the tables explicitly, and a null argument means + // "fall back to the global" so un-migrated callers keep working. Query paths that pass + // their own tables need no set_active_breakpoints() call at all. extern const float *daisy_active_breakpoints; extern const float *daisy_active_breakpoints_max; @@ -34,7 +40,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp = nullptr); float minidist_paa_to_isax_raw_SIMD(float *paa, sax_type *sax, sax_type *sax_cardinalities, @@ -43,7 +50,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp = nullptr); float ts_euclidean_distance(ts_type *t, ts_type *s, int size, float bound); @@ -56,7 +64,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp = nullptr); float minidist_paa_to_isax_raw_DTW_SING_SIMD(float *paaU, float *paaL, sax_type *sax, sax_type *sax_cardinalities, @@ -65,7 +74,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp_max = nullptr); float minidist_paa_to_isax_raw_DTW_SIMD(float *paaU, float *paaL, sax_type *sax, sax_type *sax_cardinalities, @@ -74,7 +84,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp = nullptr); float lb_keogh_data_bound(float *qo, float *tu, float *tl, float *cb, int len, float bsf); @@ -87,7 +98,8 @@ namespace daisy int number_of_segments, int min_val, int max_val, - float ratio_sqrt); + float ratio_sqrt, + const float *bp = nullptr); } diff --git a/tests/test_equidepth.cpp b/tests/test_equidepth.cpp index 6d2b265..3c11a40 100644 --- a/tests/test_equidepth.cpp +++ b/tests/test_equidepth.cpp @@ -11,7 +11,9 @@ #include #include #include +#include #include +#include #include using daisy::idx_t; @@ -139,6 +141,122 @@ TEST_F(EquidepthTest, HerculesRejectsNonNormalized) EXPECT_THROW(algo.setNormalized(false), std::runtime_error); } +// Two live indices with different equi-depth breakpoints. +// +// Breakpoint tables used to be installed into a process-global before every search, so an +// 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. +class TwoIndexEquidepthTest : public ::testing::Test +{ +protected: + static constexpr int N = 3000; + static constexpr int DIM = 64; + static constexpr int NQ = 20; + static constexpr int K = 10; + + std::vector dataA, queryA, dataB; + + void SetUp() override + { + dataA = genNonNormalized(N, DIM, 4242); + queryA = genNonNormalized(NQ, DIM, 777); + // Shift B far away so its equi-depth breakpoints cannot stand in for A's. + dataB = genNonNormalized(N, DIM, 8888); + for (float &v : dataB) + v = v * 25.0f + 100000.0f; + } + + template + void expectSecondIndexDoesNotPerturbFirst() + { + Algo a(daisy::DistanceType::L2_SQUARED); + a.setNormalized(false); + a.buildIndex(dataA.data(), N, DIM); + a.setNumThreads(1); + + std::vector alone_I((size_t)NQ * K); + std::vector alone_D((size_t)NQ * K); + a.searchIndex(queryA.data(), NQ, K, alone_I.data(), alone_D.data()); + + // Building B installs B's tables into the legacy global. + Algo b(daisy::DistanceType::L2_SQUARED); + b.setNormalized(false); + b.buildIndex(dataB.data(), N, DIM); + b.setNumThreads(1); + + std::vector after_I((size_t)NQ * K); + std::vector after_D((size_t)NQ * K); + a.searchIndex(queryA.data(), NQ, K, after_I.data(), after_D.data()); + + for (size_t i = 0; i < alone_D.size(); i++) + { + EXPECT_EQ(after_I[i], alone_I[i]) + << "result " << i << " changed once a second index existed"; + EXPECT_FLOAT_EQ(after_D[i], alone_D[i]) + << "result " << i << " changed once a second index existed"; + } + } +}; + +// 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, LbBruteforceQueriesAreUnaffectedByASecondIndex) +{ + expectSecondIndexDoesNotPerturbFirst(); +} + +// Concurrent queries against two indices with different breakpoints. A process-global +// could not express this at all: the two searches would race to install their tables. +TEST_F(TwoIndexEquidepthTest, ConcurrentQueriesOnTwoIndicesStayStable) +{ + daisy::LbBruteforce a(daisy::DistanceType::L2_SQUARED); + a.setNormalized(false); + a.buildIndex(dataA.data(), N, DIM); + a.setNumThreads(1); + + std::vector want_I((size_t)NQ * K); + std::vector want_D((size_t)NQ * K); + a.searchIndex(queryA.data(), NQ, K, want_I.data(), want_D.data()); + + daisy::LbBruteforce b(daisy::DistanceType::L2_SQUARED); + b.setNormalized(false); + b.buildIndex(dataB.data(), N, DIM); + b.setNumThreads(1); + + constexpr int ROUNDS = 40; + std::atomic stop{false}; + + // Keep B querying so its searches overlap A's for the whole run. + std::thread noise([&] + { + std::vector I((size_t)NQ * K); + std::vector D((size_t)NQ * K); + while (!stop) + b.searchIndex(dataB.data(), NQ, K, I.data(), D.data()); + }); + + std::vector I((size_t)NQ * K); + std::vector D((size_t)NQ * K); + for (int r = 0; r < ROUNDS; r++) + { + a.searchIndex(queryA.data(), NQ, K, I.data(), D.data()); + for (size_t i = 0; i < D.size(); i++) + { + ASSERT_EQ(I[i], want_I[i]) << "round " << r << ", result " << i; + ASSERT_FLOAT_EQ(D[i], want_D[i]) << "round " << r << ", result " << i; + } + } + stop = true; + noise.join(); +} + int main(int argc, char **argv) { ::testing::InitGoogleTest(&argc, argv);