diff --git a/CMakeLists.txt b/CMakeLists.txt index e56aae1..8f75903 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -40,6 +40,7 @@ option(BUILD_ODYSSEY "Build Odyssey components with MPI support" OFF) option(ODYSSEY_DEBUG_RESULTS "Print Odyssey search results (query_result) for debug before copy to I,D" OFF) option(BUILD_SING "Enable CUDA support" OFF) option(BUILD_SOFA "Build SOFA with FFTW3 support" OFF) +option(BUILD_COCONUT "Build the COCONUT sortable-SAX algorithm (static + streaming)" ON) # Backward compatibility: ENABLE_MPI and ODYSSEY_MPI map to BUILD_ODYSSEY if(DEFINED ENABLE_MPI) @@ -74,6 +75,7 @@ message(STATUS " BUILD_TESTS: ${BUILD_TESTS}") message(STATUS " BUILD_ODYSSEY: ${BUILD_ODYSSEY}") message(STATUS " ODYSSEY_DEBUG_RESULTS: ${ODYSSEY_DEBUG_RESULTS}") message(STATUS " BUILD_SING: ${BUILD_SING}") +message(STATUS " BUILD_COCONUT: ${BUILD_COCONUT}") # Define PROJECT_ROOT_DIR for all targets (paramSetup.cpp is compiled by commons_lib and by benchmark executables) get_filename_component(PROJECT_ROOT_DIR "${CMAKE_SOURCE_DIR}" ABSOLUTE) diff --git a/README.md b/README.md index b7bce4c..ff727c2 100644 --- a/README.md +++ b/README.md @@ -48,6 +48,7 @@ The following table summarizes the key features of each algorithm: | **[Hercules](https://helios2.mi.parisdescartes.fr/~themisp/publications/pvldb22-hercules.pdf)** | In-memory hierarchical similarity search using EAPCA and SAX-based pruning | | **[DumpyOS](https://helios2.mi.parisdescartes.fr/~themisp/publications/vldbj24-dumpyos.pdf)** | In-memory scalable data series similarity search using an adaptive multi-ary iSAX index | | **[FreSH](http://publications.ics.forth.gr/tech-reports/2023/2023.TR489_FreSh_A_LockFree_Data_Series_Index.pdf)** | In-memory lock-free parallel similarity search using an iSAX index (SRDS 2023) | +| **[COCONUT](http://www.vldb.org/pvldb/vol11/p677-kondylakis.pdf)** | Sortable-SAX index built bottom-up; supports both static datasets and **streaming** (incremental) inserts (PVLDB 2018) | @@ -83,6 +84,7 @@ Based on the available hardware, you can specify the below arguments to enable/d | `BUILD_DEMO` | Build demonstration applications | `ON` | Core library | | `BUILD_ODYSSEY` | Enable MPI for distributed computing | `OFF` | OpenMPI/MPICH | | `BUILD_SING` | Enable CUDA for GPU acceleration | `OFF` | CUDA Toolkit | +| `BUILD_COCONUT` | Enable the COCONUT sortable-SAX algorithm | `ON` | Core library | | `DEBUG_MSG` | Enable debug output | `OFF` | None | | `BUILD_SOFA` | Enable SOFA with SFA-Based indexing | `OFF` | FFTW3 | @@ -197,4 +199,3 @@ The logo of DaiSy was designed by [Eva Chamilaki](https://evachamilaki.github.io - diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt index 11bffd0..05f951f 100644 --- a/benchmark/CMakeLists.txt +++ b/benchmark/CMakeLists.txt @@ -140,6 +140,30 @@ if(DEBUG_MSG) message(STATUS "Include directories added for bm_LbBruteforce_L2Square.") endif() +# ////// COCONUT ////// +if(BUILD_COCONUT) + add_executable(bm_Coconut_L2Square + bm_Coconut_L2Square.cpp + bm_utils.cpp + ../commons/paramSetup.cpp + ../commons/test_bm_utils.cpp + ../commons/dataloaders.cpp + ) + target_link_libraries(bm_Coconut_L2Square + PRIVATE + benchmark::benchmark + benchmark::benchmark_main + dino_lib + ) + target_include_directories(bm_Coconut_L2Square + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) +elseif(DEBUG_MSG) + message(STATUS "BUILD_COCONUT is OFF. Skipping Coconut benchmarks.") +endif() + # ////// MESSI ////// if(DEBUG_MSG) message(STATUS "---") @@ -590,4 +614,4 @@ endif() if(DEBUG_MSG) message(STATUS "--- Benchmark Configuration Complete ---") -endif() \ No newline at end of file +endif() diff --git a/benchmark/bm_Coconut_L2Square.cpp b/benchmark/bm_Coconut_L2Square.cpp new file mode 100644 index 0000000..315246a --- /dev/null +++ b/benchmark/bm_Coconut_L2Square.cpp @@ -0,0 +1,18 @@ +#include +#include "bm_utils.hpp" +#include "../lib/algos/Coconut.hpp" + +static void BM_Coconut(benchmark::State& state) { + int config_idx = static_cast(state.range(0)); + const SSTestConfig& config = test_configs[config_idx]; + + daisy::Coconut search(daisy::DistanceType::L2_SQUARED); + + for (auto _ : state) { + runSSTBenchmark(&search, config.dataset_path, config.query_path, config.thread_count, config.k_value); + } +} + +BENCHMARK(BM_Coconut)->Arg(0)->MinTime(2.0)->Unit(benchmark::kMillisecond); + +BENCHMARK_MAIN(); diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt index 2d500ec..f4cc45f 100644 --- a/demos/CMakeLists.txt +++ b/demos/CMakeLists.txt @@ -37,6 +37,27 @@ if(BUILD_DEMO) message(STATUS "Include directories added for demo_bruteforce_L2Square.") endif() + # ////// COCONUT (static + streaming) ////// + if(BUILD_COCONUT) + if(DEBUG_MSG) + message(STATUS "## Demo: Coconut") + endif() + add_executable(demo_Coconut_L2Square demo_Coconut_L2Square.cpp) + target_link_libraries(demo_Coconut_L2Square PRIVATE dino_lib commons_lib) + target_include_directories(demo_Coconut_L2Square PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + add_executable(demo_Coconut_Streaming demo_Coconut_Streaming.cpp) + target_link_libraries(demo_Coconut_Streaming PRIVATE dino_lib commons_lib) + target_include_directories(demo_Coconut_Streaming PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + elseif(DEBUG_MSG) + message(STATUS "BUILD_COCONUT is OFF. Skipping Coconut demos.") + endif() + # ////// BRUTEFORCE DTW ////// if(DEBUG_MSG) message(STATUS "---") @@ -562,4 +583,3 @@ else() message(STATUS "BUILD_DEMO is FALSE. Demo executables will NOT be built.") endif() endif() - diff --git a/demos/demo_Coconut_L2Square.cpp b/demos/demo_Coconut_L2Square.cpp new file mode 100644 index 0000000..e8a6e5f --- /dev/null +++ b/demos/demo_Coconut_L2Square.cpp @@ -0,0 +1,43 @@ +// COCONUT — static (batch) usage: build the sortable-SAX index bottom-up from a whole +// dataset, then run exact kNN with L2 (squared) distance. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +int main() +{ + daisy::idx_t n_database = 200000; + unsigned long long dim = 96; + unsigned long long n_query = 10; + daisy::idx_t k = 5; + + 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::Coconut coconut(daisy::DistanceType::L2_SQUARED); + + // Static bottom-up build over the full dataset. + coconut.buildIndex(database, n_database, dim); + + daisy::idx_t *I = new daisy::idx_t[n_query * k]; + float *D = new float[n_query * k]; + coconut.searchIndex(query, n_query, k, I, D); + + for (daisy::idx_t i = 0; i < n_query; i++) + { + printf("Query %llu: ", i); + for (daisy::idx_t j = 0; j < k; j++) + printf("%llu ", I[i * k + j]); + printf("\n"); + } + + delete[] database; + delete[] query; + delete[] I; + delete[] D; + + return 0; +} diff --git a/demos/demo_Coconut_L2Square.py b/demos/demo_Coconut_L2Square.py new file mode 100644 index 0000000..548babe --- /dev/null +++ b/demos/demo_Coconut_L2Square.py @@ -0,0 +1,34 @@ +import sys +import os +import numpy as np + +sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) + +from daisy import DistanceType, Coconut + +def main(): + n_database = 200000 + dim = 96 + n_query = 10 + k = 5 + + np.random.seed(100) + db = np.random.randn(n_database, dim).astype(np.float32) + + np.random.seed(50) + query = np.random.randn(n_query, dim).astype(np.float32) + + index = Coconut(DistanceType.L2_SQUARED) + index.setNumThreads(1) + index.buildIndex(db) # static bottom-up build + + I, D = index.searchIndex(query, k) + + for query_num in range(n_query): + print(f"Query {query_num}:") + print("Distances:", D[query_num]) + print("Indices:", I[query_num]) + print() + +if __name__ == "__main__": + main() diff --git a/demos/demo_Coconut_Streaming.cpp b/demos/demo_Coconut_Streaming.cpp new file mode 100644 index 0000000..9ae059b --- /dev/null +++ b/demos/demo_Coconut_Streaming.cpp @@ -0,0 +1,56 @@ +// COCONUT — streaming: build on an initial batch, then insert new series into the live +// index as they arrive and query at any point (no rebuild). COCONUT's distinguishing feature. + +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" + +int main() +{ + unsigned long long dim = 96; + daisy::idx_t initial = 50000; // series available at build time + daisy::idx_t batch = 25000; // series that arrive later, per streaming step + int steps = 3; + unsigned long long n_query = 5; + daisy::idx_t k = 5; + + // Whole stream in one buffer; we feed it in pieces to simulate arrival over time. + daisy::idx_t total = initial + (daisy::idx_t)steps * batch; + float *stream = loadRandomData(total, dim, 100, true); + float *query = loadRandomData(n_query, dim, 50, true); + + daisy::Coconut coconut(daisy::DistanceType::L2_SQUARED); + + // 1. Build on the initial batch. + coconut.buildIndex(stream, initial, dim); + printf("Built on %llu series.\n", initial); + + daisy::idx_t *I = new daisy::idx_t[n_query * k]; + float *D = new float[n_query * k]; + + auto run_query = [&](const char *when) + { + coconut.searchIndex(query, n_query, k, I, D); + printf("[%s] query 0 kNN: ", when); + for (daisy::idx_t j = 0; j < k; j++) + printf("%llu(%.1f) ", I[j], D[j]); + printf("\n"); + }; + run_query("after build"); + + // 2. Stream the remaining series in batches, querying after each. + daisy::idx_t seen = initial; + for (int s = 0; s < steps; s++) + { + coconut.insertBatch(stream + (size_t)seen * dim, batch); + seen += batch; + printf("Streamed a batch; index now holds %llu series.\n", seen); + run_query("after insert"); + } + + delete[] stream; + delete[] query; + delete[] I; + delete[] D; + + return 0; +} diff --git a/demos/demo_Coconut_Streaming.py b/demos/demo_Coconut_Streaming.py new file mode 100644 index 0000000..b684267 --- /dev/null +++ b/demos/demo_Coconut_Streaming.py @@ -0,0 +1,43 @@ +import sys +import os +import numpy as np + +sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) + +from daisy import DistanceType, Coconut + +# COCONUT streaming: build on an initial batch, then insert new series into the live index +# as they arrive and query after each step (no rebuild). +def main(): + dim = 96 + initial = 50000 + batch = 25000 + steps = 3 + n_query = 5 + k = 5 + + total = initial + steps * batch + np.random.seed(100) + stream = np.random.randn(total, dim).astype(np.float32) + np.random.seed(50) + query = np.random.randn(n_query, dim).astype(np.float32) + + index = Coconut(DistanceType.L2_SQUARED) + index.setNumThreads(1) + index.buildIndex(stream[:initial]) # build on the initial batch + + def show(when): + I, D = index.searchIndex(query, k) + print(f"[{when}] query 0 -> indices {I[0]}, distances {D[0]}") + + show("after build") + + seen = initial + for s in range(steps): + index.insertBatch(stream[seen:seen + batch]) # stream the next batch + seen += batch + print(f"streamed a batch; index now holds {seen} series") + show("after insert") + +if __name__ == "__main__": + main() diff --git a/docs/demos-guide.md b/docs/demos-guide.md index 637c11a..32b3315 100644 --- a/docs/demos-guide.md +++ b/docs/demos-guide.md @@ -3,6 +3,11 @@ The demos module provides practical examples of how to use the DaiSy library's algorithms for data series similarity search. 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(...)`. +The **Coconut** algorithm additionally supports **streaming**: `demo_Coconut_L2Square` shows +the static (batch) build, while `demo_Coconut_Streaming` builds on an initial batch and then +`insert`/`insertBatch`es new series into the live index, querying after each step. + ## Demo Program Structure All demos are located in the [`demos/`](../demos/) directory. diff --git a/docs/how-to-contribute.md b/docs/how-to-contribute.md index 6afd83e..0182b02 100644 --- a/docs/how-to-contribute.md +++ b/docs/how-to-contribute.md @@ -11,7 +11,8 @@ You may reach us by contacting the main authors of the project, [Francesca Del G 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, Coconut, Dumpy, etc.) +- Extension of DaiSy for more algorithms (e.g., SFA, Hercules, Dumpy, etc.) +- Streaming / updatable indexing for more algorithms (currently supported by 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/CMakeLists.txt b/lib/algos/CMakeLists.txt index 2981ef5..890bfb4 100644 --- a/lib/algos/CMakeLists.txt +++ b/lib/algos/CMakeLists.txt @@ -27,6 +27,22 @@ if(DEBUG_MSG) message(STATUS "Bruteforce.cpp added.") endif() +# ////// COCONUT ////// +if(BUILD_COCONUT) + if(DEBUG_MSG) + message(STATUS "Adding Coconut.cpp to dino_lib sources.") + endif() + target_sources(dino_lib + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/Coconut.cpp + ) + target_compile_definitions(dino_lib PUBLIC COCONUT_ENABLED) +else() + if(DEBUG_MSG) + message(STATUS "BUILD_COCONUT is OFF. Skipping Coconut sources.") + endif() +endif() + # ////// LBBRUTEFORCE ////// if(DEBUG_MSG) message(STATUS "Adding LbBruteforce.cpp to dino_lib sources.") diff --git a/lib/algos/Coconut.cpp b/lib/algos/Coconut.cpp new file mode 100644 index 0000000..44a57d0 --- /dev/null +++ b/lib/algos/Coconut.cpp @@ -0,0 +1,266 @@ +#include "Coconut.hpp" + +#include "../isax/SAX.hpp" +#include "../isax/iSAXTypes.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace daisy +{ + + Coconut::Coconut(DistanceType distance_type) + : SimilaritySearchAlgorithm(distance_type) + { + } + + Coconut::~Coconut() + { + cleanupLeafFiles(); + } + + std::string Coconut::leafPath(int leaf_no) const + { + return index_dir_ + "/leaf." + std::to_string(leaf_no); + } + + // Compare two sortable-SAX keys lexicographically over their bytes (COCONUT's order). + static bool key_less(const std::vector &a, const std::vector &b) + { + for (size_t i = 0; i < a.size(); i++) + { + if (a[i] != b[i]) + return a[i] < b[i]; + } + return false; + } + + void Coconut::buildIndex(DataSource *data_source) + { + this->dim = data_source->getDim(); + this->n_database = data_source->getTotalRecords(); + if (this->n_database == 0) + { + data_source->reset(); + idx_t count = 0; + float *dummy = new float[this->dim]; + while (data_source->nextRecord(dummy)) + count++; + delete[] dummy; + this->n_database = count; + data_source->reset(); + } + + const int seg = this->paa_segments; + const int bits = this->sax_cardinality; + this->ts_values_per_segment_ = (int)this->dim / seg; + this->sax_alphabet_cardinality_ = (int)std::pow(2, bits); + this->mindist_sqrt_ = (float)this->dim / (float)seg; + this->max_cardinalities_.assign(seg, (sax_type)bits); + + // Read the whole database into a temporary buffer (freed once leaves are on disk). + std::vector buffer((size_t)this->n_database * this->dim); + { + float *rec = new float[this->dim]; + idx_t idx = 0; + while (data_source->nextRecord(rec) && idx < this->n_database) + { + std::copy(rec, rec + this->dim, buffer.data() + (size_t)idx * this->dim); + idx++; + } + delete[] rec; + } + + // 1. Compute SAX + sortable SAX for every series. + records_.resize(this->n_database); + for (idx_t i = 0; i < this->n_database; i++) + { + Record &r = records_[i]; + r.series_id = i; + r.sax.resize(seg); + r.key.resize(seg); + sax_from_ts(buffer.data() + (size_t)i * this->dim, r.sax.data(), + this->ts_values_per_segment_, seg, + this->sax_alphabet_cardinality_, bits); + sortable_sax_from_sax(r.sax.data(), r.key.data(), seg, bits); + } + + // 2. Sort records by sortable-SAX key (bottom-up ordering). + std::sort(records_.begin(), records_.end(), + [](const Record &a, const Record &b) { return key_less(a.key, b.key); }); + + // 3. Write on-disk leaf files: consecutive sorted records, each = [inv-SAX][raw TS]. + static std::atomic uid{0}; + this->index_dir_ = "/tmp/daisy_coconut_" + std::to_string((unsigned long)getpid()) + + "_" + std::to_string(uid.fetch_add(1)); + mkdir(this->index_dir_.c_str(), 0777); + + FILE *lf = nullptr; + int cur_leaf = -1; + for (size_t p = 0; p < records_.size(); p++) + { + int leaf_no = (int)(p / (size_t)this->leaf_capacity); + if (leaf_no != cur_leaf) + { + if (lf) fclose(lf); + lf = fopen(leafPath(leaf_no).c_str(), "wb"); + cur_leaf = leaf_no; + } + if (lf) + { + fwrite(records_[p].key.data(), sizeof(sax_type), seg, lf); + fwrite(buffer.data() + (size_t)records_[p].series_id * this->dim, + sizeof(float), this->dim, lf); + } + } + if (lf) fclose(lf); + // buffer freed here: raw series now live on disk in the leaf files. + + this->next_id_ = this->n_database; // streaming inserts get ids after the built set + } + + void Coconut::insert(const float *series) + { + if (this->paa_segments <= 0 || this->ts_values_per_segment_ <= 0) + throw std::runtime_error("Coconut::insert requires an initial buildIndex first"); + + const int seg = this->paa_segments; + Record r; + r.series_id = this->next_id_++; + r.sax.resize(seg); + r.key.resize(seg); + sax_from_ts(const_cast(series), r.sax.data(), + this->ts_values_per_segment_, seg, + this->sax_alphabet_cardinality_, this->sax_cardinality); + sortable_sax_from_sax(r.sax.data(), r.key.data(), seg, this->sax_cardinality); + + buffer_records_.push_back(std::move(r)); + buffer_data_.insert(buffer_data_.end(), series, series + this->dim); + this->n_database++; + } + + void Coconut::insertBatch(const float *data, idx_t n) + { + buffer_data_.reserve(buffer_data_.size() + (size_t)n * this->dim); + for (idx_t i = 0; i < n; i++) + insert(data + (size_t)i * this->dim); + } + + void Coconut::readSeriesFromLeaf(int leaf_no, int slot, float *out) const + { + FILE *lf = fopen(leafPath(leaf_no).c_str(), "rb"); + if (!lf) + return; + const long record_bytes = (long)this->paa_segments * (long)sizeof(sax_type) + + (long)this->dim * (long)sizeof(float); + // seek past this slot's inv-SAX to its raw TS + fseek(lf, (long)slot * record_bytes + (long)this->paa_segments * (long)sizeof(sax_type), SEEK_SET); + size_t got = fread(out, sizeof(float), this->dim, lf); + (void)got; + fclose(lf); + } + + void Coconut::searchIndex(const float *query, const idx_t n_query, const idx_t k, idx_t *I, float *D) + { + const int seg = this->paa_segments; + std::vector paa(seg); + std::vector ts(this->dim); + + for (idx_t q = 0; q < n_query; q++) + { + const float *qts = query + (size_t)q * this->dim; + paa_from_ts(qts, paa.data(), seg, this->ts_values_per_segment_); + + // k-BSF as a max-heap on distance (top = current worst of the best k). + std::priority_queue> bsf; + + for (size_t p = 0; p < records_.size(); p++) + { + float kth = (bsf.size() >= (size_t)k) ? bsf.top().first : FLT_MAX; + + float mindist = minidist_paa_to_isax( + paa.data(), + const_cast(records_[p].sax.data()), + const_cast(this->max_cardinalities_.data()), + (sax_type)this->sax_cardinality, + this->sax_alphabet_cardinality_, + seg, MINVAL, MAXVAL, this->mindist_sqrt_); + + if (mindist >= kth) + continue; // lower bound already worse than the k-th best -> prune + + int leaf_no = (int)(p / (size_t)this->leaf_capacity); + int slot = (int)(p % (size_t)this->leaf_capacity); + readSeriesFromLeaf(leaf_no, slot, ts.data()); + + float dist = ts_euclidean_distance_SIMD(const_cast(qts), ts.data(), + (int)this->dim, kth); + if (bsf.size() < (size_t)k) + bsf.push({dist, records_[p].series_id}); + else if (dist < bsf.top().first) + { + bsf.pop(); + bsf.push({dist, records_[p].series_id}); + } + } + + // Scan the in-memory streaming buffer (raw TS already in memory). + for (size_t b = 0; b < buffer_records_.size(); b++) + { + float kth = (bsf.size() >= (size_t)k) ? bsf.top().first : FLT_MAX; + float mindist = minidist_paa_to_isax( + paa.data(), + const_cast(buffer_records_[b].sax.data()), + const_cast(this->max_cardinalities_.data()), + (sax_type)this->sax_cardinality, + this->sax_alphabet_cardinality_, + seg, MINVAL, MAXVAL, this->mindist_sqrt_); + if (mindist >= kth) + continue; + float *bts = buffer_data_.data() + b * (size_t)this->dim; + float dist = ts_euclidean_distance_SIMD(const_cast(qts), bts, + (int)this->dim, kth); + if (bsf.size() < (size_t)k) + bsf.push({dist, buffer_records_[b].series_id}); + else if (dist < bsf.top().first) + { + bsf.pop(); + bsf.push({dist, buffer_records_[b].series_id}); + } + } + + // Emit in ascending distance order. + idx_t got = (idx_t)bsf.size(); + for (idx_t j = got; j > 0; j--) + { + D[q * k + (j - 1)] = bsf.top().first; + I[q * k + (j - 1)] = bsf.top().second; + bsf.pop(); + } + for (idx_t j = got; j < k; j++) // pad if fewer than k available + { + D[q * k + j] = FLT_MAX; + I[q * k + j] = 0; + } + } + } + + void Coconut::cleanupLeafFiles() + { + if (index_dir_.empty()) + return; + int nleaves = (int)((records_.size() + leaf_capacity - 1) / (size_t)std::max(1, leaf_capacity)); + for (int l = 0; l < nleaves; l++) + remove(leafPath(l).c_str()); + rmdir(index_dir_.c_str()); + index_dir_.clear(); + } + +} diff --git a/lib/algos/Coconut.hpp b/lib/algos/Coconut.hpp new file mode 100644 index 0000000..7e23304 --- /dev/null +++ b/lib/algos/Coconut.hpp @@ -0,0 +1,81 @@ +#ifndef COCONUT_HPP +#define COCONUT_HPP + +#include "SimilaritySearchAlgorithm.hpp" + +#include "../isax/iSAXTypes.hpp" + +#include +#include +#include + +namespace daisy +{ + + // COCONUT (Kondylakis et al., VLDB'18): a "sortable SAX" key (bit-interleaved SAX) orders + // series so the index can be built bottom-up by sorting and extended by streaming inserts. + // buildIndex sorts records onto on-disk leaves ([inv-SAX][raw TS] each); searchIndex does + // exact kNN via SAX MINDIST + L2 refinement. Follow-ups: paged B-tree, out-of-core sort, + // LSM merge of the insert buffer, equi-depth breakpoints. + class Coconut : public SimilaritySearchAlgorithm + { + public: + Coconut(DistanceType distance_type); + + void setNumThreads(int num_threads) override { this->num_threads = num_threads; } + int getNumThreads() const { return num_threads; } + + int getPaaSegments() const { return paa_segments; } + void setPaaSegments(int v) { paa_segments = v; } + int getSaxCardinality() const { return sax_cardinality; } + void setSaxCardinality(int v) { sax_cardinality = v; } + int getLeafSize() const { return leaf_capacity; } + void setLeafSize(int v) { leaf_capacity = v; } + + using SimilaritySearchAlgorithm::buildIndex; + void buildIndex(DataSource *data_source) override; + void buildIndex(const std::string &filename, idx_t dim, idx_t n_database = 0) override + { + throw std::runtime_error("Coconut requires in-memory data. Use buildIndex(database, n_database, dim)."); + } + + void searchIndex(const float *query, const idx_t n_query, const idx_t k, idx_t *I, float *D) override; + + // Streaming: add series to a live index (needs buildIndex first). New series go into + // an in-memory buffer that queries scan alongside the on-disk index. + void insert(const float *series) override; + void insertBatch(const float *data, idx_t n) override; + + ~Coconut() override; + + private: + // After sorting, record p lives in leaf p/leaf_capacity at slot p%leaf_capacity. + struct Record + { + std::vector sax; // SAX word (for MINDIST) + std::vector key; // sortable SAX (defines the order) + idx_t series_id; // original index in the input + }; + + std::vector records_; // main index, sorted by sortable-SAX key (TS on disk) + std::string index_dir_; // unique per-instance dir for leaf files + int leaf_capacity = 1000; // records per leaf file + + // Streaming buffer: inserted series kept in memory (summaries + raw TS), scanned by queries. + std::vector buffer_records_; + std::vector buffer_data_; // buffer_records_.size() * dim + idx_t next_id_ = 0; + + int ts_values_per_segment_ = 0; // dim / paa_segments + int sax_alphabet_cardinality_ = 0; // 2^sax_cardinality + float mindist_sqrt_ = 0.0f; // (dim / paa_segments), the MINDIST scale factor + std::vector max_cardinalities_; // all = sax_cardinality + + std::string leafPath(int leaf_no) const; + void readSeriesFromLeaf(int leaf_no, int slot, float *out) const; + void cleanupLeafFiles(); + }; + +} + +#endif diff --git a/lib/algos/SimilaritySearchAlgorithm.hpp b/lib/algos/SimilaritySearchAlgorithm.hpp index 2088b45..a41d051 100644 --- a/lib/algos/SimilaritySearchAlgorithm.hpp +++ b/lib/algos/SimilaritySearchAlgorithm.hpp @@ -103,6 +103,21 @@ namespace daisy buildIndex(&data_source); } + // Streaming API: incrementally add series to a live index. Only some algorithms + // support it (currently Coconut); the default throws. + virtual void insert(const float *series) + { + (void)series; + throw std::runtime_error("streaming insert not supported by this algorithm"); + } + + virtual void insertBatch(const float *data, idx_t n) + { + (void)data; + (void)n; + throw std::runtime_error("streaming insert not supported by this algorithm"); + } + protected: bool validateSearchParams(const idx_t k, const idx_t n_query) const { diff --git a/lib/daisy.hpp b/lib/daisy.hpp index b4da494..4c90749 100644 --- a/lib/daisy.hpp +++ b/lib/daisy.hpp @@ -3,6 +3,7 @@ #include "algos/Bruteforce.hpp" #include "algos/LbBruteforce.hpp" +#include "algos/Coconut.hpp" #include "algos/Messi.hpp" #include "algos/ParIS.hpp" #include "algos/Sing.hpp" diff --git a/lib/isax/SAX.cpp b/lib/isax/SAX.cpp index 5a6b69d..c5286f0 100644 --- a/lib/isax/SAX.cpp +++ b/lib/isax/SAX.cpp @@ -7,6 +7,30 @@ namespace daisy { + // COCONUT sortable SAX (invSax): bit-interleave the SAX across segments, MSB plane first, + // so byte order preserves SAX neighbourhood. sax/out are `segments` bytes each. + void sortable_sax_from_sax(const sax_type *sax, sax_type *out, int segments, int bit_cardinality) + { + for (int j = 0; j < segments; j++) + out[j] = 0; + + int segi = 0; + int invj = bit_cardinality - 1; + for (int i = bit_cardinality - 1; i >= 0; i--) + { + for (int j = 0; j < segments; j++) + { + unsigned int bit = ((unsigned int)sax[j] >> i) & 1u; + out[segi] |= (sax_type)(bit << invj); + if (--invj < 0) + { + segi++; + invj = bit_cardinality - 1; + } + } + } + } + // Default to the Gaussian tables; replaced per-index via set_active_breakpoints(). const float *daisy_active_breakpoints = sax_breakpoints; const float *daisy_active_breakpoints_max = sax_breakpointsnew3; diff --git a/lib/isax/SAX.hpp b/lib/isax/SAX.hpp index f8c7405..7af0984 100644 --- a/lib/isax/SAX.hpp +++ b/lib/isax/SAX.hpp @@ -23,6 +23,10 @@ namespace daisy enum response paa_from_ts(const ts_type *ts_in, ts_type *paa_out, int segments, int ts_values_per_segment); enum response sax_from_paa(ts_type *paa, sax_type *sax, int segments, int cardinality, int bit_cardinality); + // COCONUT "sortable SAX": bit-interleave a SAX word so byte order preserves SAX + // neighbourhood (enables bottom-up build + ordered inserts). out/sax are `segments` bytes. + void sortable_sax_from_sax(const sax_type *sax, sax_type *out, int segments, int bit_cardinality); + float minidist_paa_to_isax(float *paa, sax_type *sax, sax_type *sax_cardinalities, sax_type max_bit_cardinality, diff --git a/pybinds/setup.cpp b/pybinds/setup.cpp index 934eac7..7a4193c 100644 --- a/pybinds/setup.cpp +++ b/pybinds/setup.cpp @@ -12,6 +12,7 @@ #include "../lib/distance_computers/DistanceComputer.hpp" #include "../lib/algos/Bruteforce.hpp" #include "../lib/algos/LbBruteforce.hpp" +#include "../lib/algos/Coconut.hpp" #include "../lib/algos/Messi.hpp" #if ODYSSEY_MPI #include "../lib/algos/hodyssey/Odyssey.hpp" @@ -257,6 +258,67 @@ PYBIND11_MODULE(_core, m) self.searchIndex(static_cast(query_buf.ptr), n_query, config, I, D); return pybind11::make_tuple(I, D); }, "Search using SearchConfig (top-k or range) and return (indices, distances)"); + ////// COCONUT ////// + pybind11::class_(m, "Coconut", "COCONUT sortable-SAX index (static bottom-up build + streaming insert)") + .def(pybind11::init(), "Create a new Coconut with the given distance metric") + + .def("getPaaSegments", &daisy::Coconut::getPaaSegments, "Get the number of PAA segments") + .def("getSaxCardinality", &daisy::Coconut::getSaxCardinality, "Get the SAX bit cardinality") + .def("getLeafSize", &daisy::Coconut::getLeafSize, "Get the records-per-leaf-file capacity") + .def("setNumThreads", &daisy::Coconut::setNumThreads, "Set the number of threads to use") + .def("setPaaSegments", &daisy::Coconut::setPaaSegments, "Set the number of PAA segments (must divide the series length)") + .def("setSaxCardinality", &daisy::Coconut::setSaxCardinality, "Set the SAX bit cardinality") + .def("setLeafSize", &daisy::Coconut::setLeafSize, "Set the records-per-leaf-file capacity") + + // Static bottom-up build from a 2D float32 numpy array. + .def("buildIndex", [](daisy::Coconut &self, pybind11::array_t db) + { + pybind11::buffer_info buf = db.request(); + if (buf.ndim != 2) + throw std::runtime_error("Database array must be 2D"); + daisy::idx_t n = buf.shape[0]; + daisy::idx_t d = buf.shape[1]; + daisy::InMemoryDataSource data_source(static_cast(buf.ptr), n, d); + self.buildIndex(&data_source); }, "Build the index from a 2D float32 numpy array") + + // Streaming: insert one series (1D) or a batch (2D) into a live index. + .def("insert", [](daisy::Coconut &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"); + self.insert(static_cast(buf.ptr)); }, "Insert a single series (1D float32 array) into the live index") + + .def("insertBatch", [](daisy::Coconut &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"); + self.insertBatch(static_cast(buf.ptr), buf.shape[0]); }, "Insert a batch of series (2D float32 array) into the live index") + + // kNN search returning (indices, distances). + .def("searchIndex", [](daisy::Coconut &self, pybind11::array_t query, daisy::idx_t k) + { + pybind11::buffer_info buf = query.request(); + if (buf.ndim != 2) + throw std::runtime_error("Query array must be 2D"); + + std::vector indices(buf.shape[0] * k); + std::vector distances(buf.shape[0] * k); + self.searchIndex(static_cast(buf.ptr), buf.shape[0], k, indices.data(), distances.data()); + + return pybind11::make_tuple( + pybind11::array_t(pybind11::buffer_info( + indices.data(), sizeof(daisy::idx_t), + pybind11::format_descriptor::format(), 2, + std::vector{static_cast(buf.shape[0]), static_cast(k)}, + std::vector{static_cast(sizeof(daisy::idx_t) * k), static_cast(sizeof(daisy::idx_t))})), + pybind11::array_t(pybind11::buffer_info( + distances.data(), sizeof(float), + pybind11::format_descriptor::format(), 2, + std::vector{static_cast(buf.shape[0]), static_cast(k)}, + std::vector{static_cast(sizeof(float) * k), static_cast(sizeof(float))}))); }, "kNN search: returns (indices, distances)"); + ////// MESSI ////// pybind11::class_(m, "Messi", "MESSI (Multi-Queue Efficient SAX Similarity Index) algorithm for time series similarity search") // Constructor diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 24cd76b..4587d98 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1005,6 +1005,24 @@ if(BUILD_ODYSSEY_AVAILABLE) WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} ENVIRONMENT "PMIX_MCA_pcompress_base_silence_warning=1") endif() +# ////// COCONUT (static kNN) ////// +if(BUILD_COCONUT) + add_executable(test_Coconut_L2Square test_Coconut_L2Square.cpp test_utils.cpp) + target_link_libraries(test_Coconut_L2Square PRIVATE GTest::gtest_main dino_lib commons_lib stdc++fs) + target_include_directories(test_Coconut_L2Square PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib ${CMAKE_CURRENT_SOURCE_DIR}/../commons) + gtest_discover_tests(test_Coconut_L2Square WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}) + + # ////// COCONUT (streaming kNN) ////// + add_executable(test_Coconut_Streaming test_Coconut_Streaming.cpp test_utils.cpp) + target_link_libraries(test_Coconut_Streaming PRIVATE GTest::gtest_main dino_lib commons_lib stdc++fs) + target_include_directories(test_Coconut_Streaming PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib ${CMAKE_CURRENT_SOURCE_DIR}/../commons) + gtest_discover_tests(test_Coconut_Streaming WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}) +elseif(DEBUG_MSG) + message(STATUS "BUILD_COCONUT is OFF. Skipping Coconut tests.") +endif() + # non-normalized data add_executable( test_equidepth diff --git a/tests/test_Coconut_L2Square.cpp b/tests/test_Coconut_L2Square.cpp new file mode 100644 index 0000000..fda674e --- /dev/null +++ b/tests/test_Coconut_L2Square.cpp @@ -0,0 +1,76 @@ +// Static COCONUT: exact kNN on a z-normalized dataset must match brute-force ground truth +// (a wrong sortable-SAX / MINDIST path would drop true neighbours). + +#include "test_utils.hpp" // gtest + algorithm headers + BruteForceSearch +#include "../lib/algos/Coconut.hpp" + +#include +#include +#include +#include + +using daisy::idx_t; + +namespace +{ + // Z-normalized random-walk series (the setting SAX/Gaussian breakpoints target). + std::vector genZNorm(int n, int dim, unsigned seed) + { + std::mt19937 rng(seed); + std::normal_distribution step(0.0f, 1.0f); + std::vector data((size_t)n * dim); + for (int i = 0; i < n; i++) + { + float *s = data.data() + (size_t)i * dim; + float acc = 0.0f; + for (int j = 0; j < dim; j++) { acc += step(rng); s[j] = acc; } + double mean = 0, sd = 0; + for (int j = 0; j < dim; j++) mean += s[j]; + mean /= dim; + for (int j = 0; j < dim; j++) sd += (s[j] - mean) * (s[j] - mean); + sd = std::sqrt(sd / dim); + if (sd < 1e-6) sd = 1.0; + for (int j = 0; j < dim; j++) s[j] = (float)((s[j] - mean) / sd); + } + return data; + } + + bool distClose(double a, double b) { return std::fabs(a - b) <= 1e-8 + 1e-2 * std::fabs(b); } +} + +TEST(CoconutTest, StaticKnnMatchesBruteForce) +{ + const int N = 2000, DIM = 64, NQ = 30, K = 10; + auto data = genZNorm(N, DIM, 7); + auto query = genZNorm(NQ, DIM, 77); + + daisy::BruteForceSearch bf(daisy::DistanceType::L2_SQUARED); + bf.buildIndex(data.data(), N, DIM); + std::vector gtI((size_t)NQ * K); + std::vector gtD((size_t)NQ * K); + bf.searchIndex(query.data(), NQ, K, gtI.data(), gtD.data()); + + daisy::Coconut cc(daisy::DistanceType::L2_SQUARED); + cc.buildIndex(data.data(), N, DIM); + cc.setNumThreads(1); + std::vector I((size_t)NQ * K); + std::vector D((size_t)NQ * K); + cc.searchIndex(query.data(), NQ, K, I.data(), D.data()); + + for (int q = 0; q < NQ; q++) + { + std::vector gd(gtD.begin() + (size_t)q * K, gtD.begin() + (size_t)(q + 1) * K); + std::vector od(D.begin() + (size_t)q * K, D.begin() + (size_t)(q + 1) * K); + std::sort(gd.begin(), gd.end()); + std::sort(od.begin(), od.end()); + for (int j = 0; j < K; j++) + EXPECT_TRUE(distClose(od[j], gd[j])) + << "query " << q << " nn " << j << ": coconut " << od[j] << " vs GT " << gd[j]; + } +} + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/tests/test_Coconut_Streaming.cpp b/tests/test_Coconut_Streaming.cpp new file mode 100644 index 0000000..376f1f8 --- /dev/null +++ b/tests/test_Coconut_Streaming.cpp @@ -0,0 +1,107 @@ +// Streaming COCONUT: build on an initial batch, then insert more batches; after each step +// exact kNN must match brute-force ground truth over the data seen so far. + +#include "test_utils.hpp" +#include "../lib/algos/Coconut.hpp" + +#include +#include +#include +#include + +using daisy::idx_t; + +namespace +{ + std::vector genZNorm(int n, int dim, unsigned seed) + { + std::mt19937 rng(seed); + std::normal_distribution step(0.0f, 1.0f); + std::vector data((size_t)n * dim); + for (int i = 0; i < n; i++) + { + float *s = data.data() + (size_t)i * dim; + float acc = 0.0f; + for (int j = 0; j < dim; j++) { acc += step(rng); s[j] = acc; } + double mean = 0, sd = 0; + for (int j = 0; j < dim; j++) mean += s[j]; + mean /= dim; + for (int j = 0; j < dim; j++) sd += (s[j] - mean) * (s[j] - mean); + sd = std::sqrt(sd / dim); + if (sd < 1e-6) sd = 1.0; + for (int j = 0; j < dim; j++) s[j] = (float)((s[j] - mean) / sd); + } + return data; + } + + bool distClose(double a, double b) { return std::fabs(a - b) <= 1e-8 + 1e-2 * std::fabs(b); } + + // exact kNN distances (sorted) of `query` over the first `n_seen` series of `data` + void bruteForceKnn(const std::vector &data, int n_seen, int dim, + const float *query, int NQ, int K, std::vector> &out) + { + daisy::BruteForceSearch bf(daisy::DistanceType::L2_SQUARED); + bf.buildIndex(const_cast(data.data()), n_seen, dim); + std::vector gtI((size_t)NQ * K); + std::vector gtD((size_t)NQ * K); + bf.searchIndex(const_cast(query), NQ, K, gtI.data(), gtD.data()); + out.assign(NQ, {}); + for (int q = 0; q < NQ; q++) + { + out[q].assign(gtD.begin() + (size_t)q * K, gtD.begin() + (size_t)(q + 1) * K); + std::sort(out[q].begin(), out[q].end()); + } + } + + void checkAgainstGT(daisy::Coconut &cc, const std::vector &data, int n_seen, int dim, + const std::vector &query, int NQ, int K, const char *phase) + { + std::vector> gt; + bruteForceKnn(data, n_seen, dim, query.data(), NQ, K, gt); + + std::vector I((size_t)NQ * K); + std::vector D((size_t)NQ * K); + cc.searchIndex(query.data(), NQ, K, I.data(), D.data()); + + for (int q = 0; q < NQ; q++) + { + std::vector od(D.begin() + (size_t)q * K, D.begin() + (size_t)(q + 1) * K); + std::sort(od.begin(), od.end()); + for (int j = 0; j < K; j++) + EXPECT_TRUE(distClose(od[j], gt[q][j])) + << phase << " query " << q << " nn " << j + << ": coconut " << od[j] << " vs GT " << gt[q][j]; + } + } +} + +TEST(CoconutTest, StreamingKnnMatchesBruteForce) +{ + const int DIM = 64, NQ = 20, K = 10; + const int INITIAL = 800; + const int BATCH = 400; + const int BATCHES = 3; // total = 800 + 3*400 = 2000 + const int N = INITIAL + BATCHES * BATCH; + + auto data = genZNorm(N, DIM, 11); + auto query = genZNorm(NQ, DIM, 111); + + daisy::Coconut cc(daisy::DistanceType::L2_SQUARED); + cc.buildIndex(data.data(), INITIAL, DIM); // static bottom-up build on the first batch + cc.setNumThreads(1); + checkAgainstGT(cc, data, INITIAL, DIM, query, NQ, K, "after build"); + + int seen = INITIAL; + for (int b = 0; b < BATCHES; b++) + { + cc.insertBatch(data.data() + (size_t)seen * DIM, BATCH); // stream the next batch + seen += BATCH; + checkAgainstGT(cc, data, seen, DIM, query, NQ, K, "after insert"); + } +} + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +}