From 7b647a9f1d071eb87aec4a06fc655900db2b16e4 Mon Sep 17 00:00:00 2001 From: sophisid Date: Fri, 24 Jul 2026 14:09:20 +0200 Subject: [PATCH] added equi-depth SAX breakpoints for all approaches. Gaussian by default (for z-norm data) --- .gitignore | 1 + lib/algos/DumpyOS.cpp | 8 +- lib/algos/Fresh.cpp | 9 +- lib/algos/LbBruteforce.cpp | 9 +- lib/algos/Messi.cpp | 9 +- lib/algos/ParIS.cpp | 28 +++- lib/algos/SimilaritySearchAlgorithm.hpp | 16 ++ lib/algos/Sing.cpp | 16 +- lib/algos/Sofa.cpp | 8 +- lib/isax/SAX.cpp | 186 +++++++++++++----------- lib/isax/SAX.hpp | 8 + lib/isax/iSAXIndex.cpp | 117 ++++++++++++++- lib/isax/iSAXIndex.hpp | 24 ++- pybinds/setup.cpp | 17 +++ 14 files changed, 340 insertions(+), 116 deletions(-) diff --git a/.gitignore b/.gitignore index 445fdfc..5c02a6e 100644 --- a/.gitignore +++ b/.gitignore @@ -37,6 +37,7 @@ *.part *.full isax_file.sax +daisy_index_nodes/ # Build directories build/ diff --git a/lib/algos/DumpyOS.cpp b/lib/algos/DumpyOS.cpp index 3026a16..efd1dda 100644 --- a/lib/algos/DumpyOS.cpp +++ b/lib/algos/DumpyOS.cpp @@ -420,13 +420,13 @@ static void get_value_range(int sax_val, int bc, double* lb, double* ub) { int offset = ((cardinality - 1) * (cardinality - 2)) / 2; if (sax_val == 0) { *lb = -std::numeric_limits::max(); - *ub = sax_breakpoints[offset]; + *ub = daisy_active_breakpoints[offset]; } else if (sax_val == cardinality - 1) { - *lb = sax_breakpoints[offset + sax_val - 1]; + *lb = daisy_active_breakpoints[offset + sax_val - 1]; *ub = std::numeric_limits::max(); } else { - *lb = sax_breakpoints[offset + sax_val - 1]; - *ub = sax_breakpoints[offset + sax_val]; + *lb = daisy_active_breakpoints[offset + sax_val - 1]; + *ub = daisy_active_breakpoints[offset + sax_val]; } } diff --git a/lib/algos/Fresh.cpp b/lib/algos/Fresh.cpp index 17bf2ba..a2afa33 100644 --- a/lib/algos/Fresh.cpp +++ b/lib/algos/Fresh.cpp @@ -899,7 +899,12 @@ void Fresh::buildIndex(DataSource *data_source) this->tight_bound, 0, 1, - 1); + 1, + this->bp_mode); + + // Equi-depth breakpoints (no-op in Gaussian mode), installed as active. + compute_equidepth_breakpoints(this->index_settings, this->database, this->n_database); + activateBreakpoints(); this->index = isax_index_init_inmemory(this->index_settings); isax_index *index = this->index; @@ -1470,6 +1475,7 @@ void Fresh::searchIndex(const float *query, idx_t n_query, const SearchConfig &c std::vector> &I, std::vector> &D) { + activateBreakpoints(); if (config.type == QueryType::TOP_K) { SimilaritySearchAlgorithm::searchIndex(query, n_query, config, I, D); @@ -1513,6 +1519,7 @@ void Fresh::searchIndex(const float *query, idx_t n_query, const SearchConfig &c void Fresh::searchIndex(const float *query, idx_t n_query, idx_t k, idx_t *I, float *D) { + activateBreakpoints(); if (this->distance_type == DistanceType::L2_SQUARED) searchIndexL2Squared(query, n_query, k, I, D); else if (this->distance_type == DistanceType::DTW) diff --git a/lib/algos/LbBruteforce.cpp b/lib/algos/LbBruteforce.cpp index feaa6bf..be317b6 100644 --- a/lib/algos/LbBruteforce.cpp +++ b/lib/algos/LbBruteforce.cpp @@ -78,7 +78,12 @@ namespace daisy this->tight_bound, 0, 1, - 1); + 1, + this->bp_mode); + + // Equi-depth breakpoints (no-op in Gaussian mode), installed as active. + compute_equidepth_breakpoints(this->index_settings, this->database, this->n_database); + activateBreakpoints(); this->index = isax_index_init_inmemory(this->index_settings); @@ -105,6 +110,7 @@ namespace daisy void LbBruteforce::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); @@ -287,6 +293,7 @@ namespace daisy std::vector> &I, std::vector> &D) { + activateBreakpoints(); if (config.type == QueryType::TOP_K) { SimilaritySearchAlgorithm::searchIndex(query, n_query, config, I, D); return; diff --git a/lib/algos/Messi.cpp b/lib/algos/Messi.cpp index 7945660..939e535 100644 --- a/lib/algos/Messi.cpp +++ b/lib/algos/Messi.cpp @@ -794,7 +794,12 @@ namespace daisy this->tight_bound, 0, 1, - 1); + 1, + this->bp_mode); + + // Equi-depth breakpoints (no-op in Gaussian mode), installed as active. + compute_equidepth_breakpoints(this->index_settings, this->database, this->n_database); + activateBreakpoints(); this->index = isax_index_init_inmemory(this->index_settings); isax_index *index = this->index; @@ -1074,6 +1079,7 @@ namespace daisy std::vector> &I, std::vector> &D) { + activateBreakpoints(); if (config.type == QueryType::TOP_K) { SimilaritySearchAlgorithm::searchIndex(query, n_query, config, I, D); @@ -1117,6 +1123,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); diff --git a/lib/algos/ParIS.cpp b/lib/algos/ParIS.cpp index e48c112..3d74655 100644 --- a/lib/algos/ParIS.cpp +++ b/lib/algos/ParIS.cpp @@ -54,10 +54,28 @@ namespace daisy this->flush_limit, this->initial_fbl_size, this->total_loaded_leaves, - this->tight_bound, - 0, - 1, - 0); + this->tight_bound, + 0, + 1, + 0, + this->bp_mode); + + // Equi-depth breakpoints. ParIS streams from a raw float32 file, so sample leading + // rows from it to derive the distribution. + if (this->bp_mode == BP_EQUIDEPTH && this->n_database > 0 && this->dim > 0) + { + size_t cap_rows = (size_t)2000000 / (size_t)std::max(1, this->paa_segments); + size_t sample_rows = std::min((size_t)this->n_database, std::max((size_t)1, cap_rows)); + std::vector sample((size_t)sample_rows * (size_t)this->dim); + FILE *sf = fopen(filename, "rb"); + if (sf) + { + size_t got = fread(sample.data(), sizeof(float), sample.size(), sf); + fclose(sf); + compute_equidepth_breakpoints(this->index_settings, sample.data(), got / (size_t)this->dim); + } + } + activateBreakpoints(); this->index = isax_index_init(this->index_settings); isax_index *index = this->index; @@ -151,6 +169,7 @@ namespace daisy void ParIS::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); @@ -1366,6 +1385,7 @@ namespace daisy std::vector> &I, std::vector> &D) { + activateBreakpoints(); if (config.type == QueryType::TOP_K) { SimilaritySearchAlgorithm::searchIndex(query, n_query, config, I, D); return; diff --git a/lib/algos/SimilaritySearchAlgorithm.hpp b/lib/algos/SimilaritySearchAlgorithm.hpp index 6d29d97..3059c9c 100644 --- a/lib/algos/SimilaritySearchAlgorithm.hpp +++ b/lib/algos/SimilaritySearchAlgorithm.hpp @@ -52,6 +52,8 @@ namespace daisy float minimum_distance = FLT_MAX; int min_checked_leaves = -1; + breakpoint_mode bp_mode = BP_GAUSSIAN; + isax_index_settings *index_settings = nullptr; isax_index *index = nullptr; sax_type **db_sax_representations = nullptr; @@ -68,6 +70,20 @@ namespace daisy idx_t getDim() const { return dim; } isax_index *getIndex() const { return index; } + // 0 = Gaussian (default), 1 = equi-depth. Set before buildIndex(). + void setBreakpointMode(int mode) { bp_mode = (mode == 1) ? BP_EQUIDEPTH : BP_GAUSSIAN; } + int getBreakpointMode() const { return (int)bp_mode; } + + protected: + // Install this index's breakpoints globally; call at build and at each search entry. + void activateBreakpoints() + { + if (index_settings) + set_active_breakpoints(index_settings->breakpoints, index_settings->breakpoints_max); + } + + public: + virtual void buildIndex(DataSource *data_source) = 0; virtual void buildIndex(float *database, idx_t n_database, idx_t dim) diff --git a/lib/algos/Sing.cpp b/lib/algos/Sing.cpp index 144fb7c..8d532f6 100644 --- a/lib/algos/Sing.cpp +++ b/lib/algos/Sing.cpp @@ -105,8 +105,8 @@ namespace daisy __m256i lower_juge_nzerov_1 = _mm256_andnot_si256(lower_juge_zerov_1, vectorsignbit); __m256 minvalv = _mm256_set1_ps((float)min_val); - __m256 lsax_breakpoints_shiftv_0 = _mm256_i32gather_ps(sax_breakpointsnew3, region_lowerv_0, 4); - __m256 lsax_breakpoints_shiftv_1 = _mm256_i32gather_ps(sax_breakpointsnew3, region_lowerv_1, 4); + __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 breakpoint_lowerv_0 = (__m256)_mm256_castsi256_ps(_mm256_or_si256( _mm256_and_si256(lower_juge_zerov_0, _mm256_castps_si256(minvalv)), _mm256_and_si256(lower_juge_nzerov_0, _mm256_castps_si256(lsax_breakpoints_shiftv_0)))); @@ -114,8 +114,8 @@ namespace daisy _mm256_and_si256(lower_juge_zerov_1, _mm256_castps_si256(minvalv)), _mm256_and_si256(lower_juge_nzerov_1, _mm256_castps_si256(lsax_breakpoints_shiftv_1)))); - __m256 usax_breakpoints_shiftv_0 = _mm256_i32gather_ps(sax_breakpointsnew3, region_upperv_0, 4); - __m256 usax_breakpoints_shiftv_1 = _mm256_i32gather_ps(sax_breakpointsnew3, region_upperv_1, 4); + __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); __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)); __m256i upper_juge_nmaxv_0 = _mm256_andnot_si256(upper_juge_maxv_0, vectorsignbit); @@ -267,7 +267,7 @@ namespace daisy } else { - breakpoint_lower = sax_breakpointsnew3[region_lower]; + breakpoint_lower = daisy_active_breakpoints_max[region_lower]; } if (region_upper == (sax_type)(max_cardinality - 1)) { @@ -275,7 +275,7 @@ namespace daisy } else { - breakpoint_upper = sax_breakpointsnew3[region_upper + 1]; + breakpoint_upper = daisy_active_breakpoints_max[region_upper + 1]; } if (breakpoint_lower > paa[i]) @@ -471,7 +471,7 @@ namespace daisy } else { - breakpoint_lower = sax_breakpointsnew3[region_lower - 1]; + breakpoint_lower = daisy_active_breakpoints_max[region_lower - 1]; } if (region_upper == (sax_type)(max_cardinality - 1)) { @@ -479,7 +479,7 @@ namespace daisy } else { - breakpoint_upper = sax_breakpointsnew3[region_upper]; + breakpoint_upper = daisy_active_breakpoints_max[region_upper]; } (void)mask; diff --git a/lib/algos/Sofa.cpp b/lib/algos/Sofa.cpp index e0bd447..d7e2e56 100644 --- a/lib/algos/Sofa.cpp +++ b/lib/algos/Sofa.cpp @@ -867,6 +867,7 @@ void Sofa::searchIndexDTW(const float *query, idx_t n_query, idx_t k, idx_t *I, void Sofa::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); else if (this->distance_type == DistanceType::DTW) @@ -879,6 +880,7 @@ void Sofa::searchIndex(const float *query, idx_t n_query, const SearchConfig &co std::vector> &I, std::vector> &D) { + activateBreakpoints(); if (config.type == QueryType::TOP_K) { SimilaritySearchAlgorithm::searchIndex(query, n_query, config, I, D); return; @@ -1047,7 +1049,11 @@ void Sofa::buildIndex(DataSource *data_source) this->leaf_size, this->min_leaf_size, this->initial_lbl_size, this->flush_limit, this->initial_fbl_size, this->total_loaded_leaves, - this->tight_bound, 0, 1, 1); + this->tight_bound, 0, 1, 1, this->bp_mode); + + // Equi-depth breakpoints (no-op in Gaussian mode), installed as active. + compute_equidepth_breakpoints(this->index_settings, this->database, this->n_database); + activateBreakpoints(); this->index = isax_index_init_inmemory(this->index_settings); isax_index *index = this->index; diff --git a/lib/isax/SAX.cpp b/lib/isax/SAX.cpp index 8eb8f91..5a6b69d 100644 --- a/lib/isax/SAX.cpp +++ b/lib/isax/SAX.cpp @@ -7,6 +7,16 @@ namespace daisy { + // 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; + + void set_active_breakpoints(const float *breakpoints, const float *breakpoints_max) + { + daisy_active_breakpoints = breakpoints ? breakpoints : sax_breakpoints; + daisy_active_breakpoints_max = breakpoints_max ? breakpoints_max : sax_breakpointsnew3; + } + int compare(const void *a, const void *b) { float *c = (float *)b - 1; @@ -57,14 +67,14 @@ namespace daisy { sax_out[si] = 0; - // First object = sax_breakpoints[offset] - // Last object = sax_breakpoints[offset + cardinality - 2] + // First object = daisy_active_breakpoints[offset] + // Last object = daisy_active_breakpoints[offset + cardinality - 2] // Size of sub-array = cardinality - 1 - float *res = (float *)bsearch(&paa[si], &sax_breakpoints[offset], cardinality - 1, + float *res = (float *)bsearch(&paa[si], &daisy_active_breakpoints[offset], cardinality - 1, sizeof(ts_type), compare); if (res != NULL) - sax_out[si] = (int)(res - &sax_breakpoints[offset]); + sax_out[si] = (int)(res - &daisy_active_breakpoints[offset]); else if (paa[si] > 0) sax_out[si] = cardinality - 1; } @@ -97,14 +107,14 @@ namespace daisy { sax[si] = 0; - // First object = sax_breakpoints[offset] - // Last object = sax_breakpoints[offset + cardinality - 2] + // First object = daisy_active_breakpoints[offset] + // Last object = daisy_active_breakpoints[offset + cardinality - 2] // Size of sub-array = cardinality - 1 - float *res = (float *)bsearch(&paa[si], &sax_breakpoints[offset], cardinality - 1, + float *res = (float *)bsearch(&paa[si], &daisy_active_breakpoints[offset], cardinality - 1, sizeof(ts_type), compare); if (res != NULL) - sax[si] = (int)(res - &sax_breakpoints[offset]); + sax[si] = (int)(res - &daisy_active_breakpoints[offset]); else if (paa[si] > 0) sax[si] = cardinality - 1; } @@ -148,7 +158,7 @@ namespace daisy } else { - breakpoint_lower = sax_breakpoints[offset + region_lower - 1]; + breakpoint_lower = daisy_active_breakpoints[offset + region_lower - 1]; } if (region_upper == max_cardinality - 1) { @@ -156,7 +166,7 @@ namespace daisy } else { - breakpoint_upper = sax_breakpoints[offset + region_upper]; + breakpoint_upper = daisy_active_breakpoints[offset + region_upper]; } if (breakpoint_lower > paa[i]) @@ -243,43 +253,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(sax_breakpoints[offset + region_lower[7] - 1], - sax_breakpoints[offset + region_lower[6] - 1], - sax_breakpoints[offset + region_lower[5] - 1], - sax_breakpoints[offset + region_lower[4] - 1], - sax_breakpoints[offset + region_lower[3] - 1], - sax_breakpoints[offset + region_lower[2] - 1], - sax_breakpoints[offset + region_lower[1] - 1], - sax_breakpoints[offset + region_lower[0] - 1]); - __m256 lsax_breakpoints_shiftv_1 = _mm256_set_ps(sax_breakpoints[offset + region_lower[15] - 1], - sax_breakpoints[offset + region_lower[14] - 1], - sax_breakpoints[offset + region_lower[13] - 1], - sax_breakpoints[offset + region_lower[12] - 1], - sax_breakpoints[offset + region_lower[11] - 1], - sax_breakpoints[offset + region_lower[10] - 1], - sax_breakpoints[offset + region_lower[9] - 1], - sax_breakpoints[offset + region_lower[8] - 1]); + __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 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(sax_breakpoints[offset + region_upper[7]], - sax_breakpoints[offset + region_upper[6]], - sax_breakpoints[offset + region_upper[5]], - sax_breakpoints[offset + region_upper[4]], - sax_breakpoints[offset + region_upper[3]], - sax_breakpoints[offset + region_upper[2]], - sax_breakpoints[offset + region_upper[1]], - sax_breakpoints[offset + region_upper[0]]); - __m256 usax_breakpoints_shiftv_1 = _mm256_set_ps(sax_breakpoints[offset + region_upper[15]], - sax_breakpoints[offset + region_upper[14]], - sax_breakpoints[offset + region_upper[13]], - sax_breakpoints[offset + region_upper[12]], - sax_breakpoints[offset + region_upper[11]], - sax_breakpoints[offset + region_upper[10]], - sax_breakpoints[offset + region_upper[9]], - sax_breakpoints[offset + region_upper[8]]); + __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]]); __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)); @@ -405,8 +415,8 @@ namespace daisy __m256 minvalv = _mm256_set1_ps((float)min_val); - __m256 lsax_breakpoints_shiftv_0 = _mm256_i32gather_ps(sax_breakpointsnew3, region_lowerv_0, 4); - __m256 lsax_breakpoints_shiftv_1 = _mm256_i32gather_ps(sax_breakpointsnew3, region_lowerv_1, 4); + __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 breakpoint_lowerv_0 = (__m256)_mm256_or_si256( _mm256_and_si256(lower_juge_zerov_0, (__m256i)minvalv), @@ -415,8 +425,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(sax_breakpointsnew3, region_upperv_0, 4); - __m256 usax_breakpoints_shiftv_1 = _mm256_i32gather_ps(sax_breakpointsnew3, region_upperv_1, 4); + __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); __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)); @@ -618,15 +628,15 @@ namespace daisy __m256 minvalv = _mm256_set1_ps(min_val); - __m256 lsax_breakpoints_shiftv_0 = _mm256_i32gather_ps(sax_breakpoints, region_lowerv_0_offset, 4); - __m256 lsax_breakpoints_shiftv_1 = _mm256_i32gather_ps(sax_breakpoints, region_lowerv_1_offset, 4); + __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 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(sax_breakpoints, region_upperv_0_offset, 4); - __m256 usax_breakpoints_shiftv_1 = _mm256_i32gather_ps(sax_breakpoints, region_upperv_1_offset, 4); + __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); __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)); @@ -762,47 +772,47 @@ namespace daisy __m256 minvalv = _mm256_set1_ps(min_val); - __m256 lsax_breakpoints_shiftv_0 = _mm256_i32gather_ps(sax_breakpoints, region_lowerv_0_offset, 4); - //__m256 lsax_breakpoints_shiftv_0= _mm256_set_ps (sax_breakpoints[region_lower[7]], - // sax_breakpoints[region_lower[6]], - // sax_breakpoints[region_lower[5]], - // sax_breakpoints[region_lower[4]], - // sax_breakpoints[region_lower[3]], - // sax_breakpoints[region_lower[2]], - // sax_breakpoints[region_lower[1]], - // sax_breakpoints[region_lower[0]]); - __m256 lsax_breakpoints_shiftv_1 = _mm256_i32gather_ps(sax_breakpoints, region_lowerv_1_offset, 4); - //__m256 lsax_breakpoints_shiftv_1= _mm256_set_ps (sax_breakpoints[region_lower[15]], - // sax_breakpoints[region_lower[14]], - // sax_breakpoints[region_lower[13]], - // sax_breakpoints[region_lower[12]], - // sax_breakpoints[region_lower[11]], - // sax_breakpoints[region_lower[10]], - // sax_breakpoints[region_lower[9]], - // sax_breakpoints[region_lower[8]]); + __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 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(sax_breakpoints, region_upperv_0_offset, 4); - //__m256 usax_breakpoints_shiftv_0= _mm256_set_ps (sax_breakpoints[region_upper[7]], - // sax_breakpoints[region_upper[6]], - // sax_breakpoints[region_upper[5]], - // sax_breakpoints[region_upper[4]], - // sax_breakpoints[region_upper[3]], - // sax_breakpoints[region_upper[2]], - // sax_breakpoints[region_upper[1]], - // sax_breakpoints[region_upper[0]]); - __m256 usax_breakpoints_shiftv_1 = _mm256_i32gather_ps(sax_breakpoints, region_upperv_1_offset, 4); - //__m256 usax_breakpoints_shiftv_1= _mm256_set_ps (sax_breakpoints[region_upper[15]], - // sax_breakpoints[region_upper[14]], - // sax_breakpoints[region_upper[13]], - // sax_breakpoints[region_upper[12]], - // sax_breakpoints[region_upper[11]], - // sax_breakpoints[region_upper[10]], - // sax_breakpoints[region_upper[9]], - // sax_breakpoints[region_upper[8]]); + __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]]); __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)); @@ -1079,7 +1089,7 @@ namespace daisy } else { - breakpoint_lower = sax_breakpoints[offset + region_lower - 1]; + breakpoint_lower = daisy_active_breakpoints[offset + region_lower - 1]; } if (region_upper == max_cardinality - 1) { @@ -1087,7 +1097,7 @@ namespace daisy } else { - breakpoint_upper = sax_breakpoints[offset + region_upper]; + breakpoint_upper = daisy_active_breakpoints[offset + region_upper]; } if (breakpoint_lower > paaU[i]) diff --git a/lib/isax/SAX.hpp b/lib/isax/SAX.hpp index fa4d6b1..f8c7405 100644 --- a/lib/isax/SAX.hpp +++ b/lib/isax/SAX.hpp @@ -9,6 +9,14 @@ 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. + extern const float *daisy_active_breakpoints; + extern const float *daisy_active_breakpoints_max; + + // NULL args restore the Gaussian defaults. + void set_active_breakpoints(const float *breakpoints, const float *breakpoints_max); + int compare(const void *a, const void *b); enum response sax_from_ts(ts_type *ts_in, sax_type *sax_out, int ts_values_per_segment, int segments, int cardinality, int bit_cardinality); diff --git a/lib/isax/iSAXIndex.cpp b/lib/isax/iSAXIndex.cpp index 0b996a5..a23a152 100644 --- a/lib/isax/iSAXIndex.cpp +++ b/lib/isax/iSAXIndex.cpp @@ -5,6 +5,9 @@ #include #include #include // for memset +#include // mkdir +#include +#include namespace daisy { @@ -13,7 +16,8 @@ namespace daisy int max_leaf_size, int min_leaf_size, int initial_leaf_buffer_size, int max_total_buffer_size, int initial_fbl_buffer_size, - int total_loaded_leaves, int tight_bound, int aggressive_check, int new_index, char inmemory_flag) + int total_loaded_leaves, int tight_bound, int aggressive_check, int new_index, char inmemory_flag, + breakpoint_mode bp_mode) { int i; isax_index_settings *settings = (isax_index_settings *)malloc(sizeof(isax_index_settings)); @@ -104,7 +108,8 @@ namespace daisy float c_size = ceil(log10(settings->sax_alphabet_cardinality + 1)); settings->max_filename_size = settings->paa_segments * ((c_size * 2) + 2) + - 5 + strlen(root_directory); + 5 + strlen(root_directory) + + 64; // margin for the default node-file directory prefix if (paa_segments > sax_bit_cardinality) { @@ -146,9 +151,95 @@ namespace daisy settings->max_total_buffer_size = settings->max_total_full_buffer_size; } + // Default to the Gaussian tables; replaced by compute_equidepth_breakpoints(). + settings->bp_mode = bp_mode; + settings->breakpoints = sax_breakpoints; + settings->breakpoints_max = sax_breakpointsnew3; + settings->breakpoints_owned = NULL; + settings->breakpoints_max_owned = NULL; + return settings; } + void compute_equidepth_breakpoints(isax_index_settings *settings, + const float *database, size_t n_series) + { + // No-op unless equi-depth was requested; keeps the Gaussian defaults. + if (settings == NULL || settings->bp_mode != BP_EQUIDEPTH) + return; + if (database == NULL || n_series == 0) + return; + + const int seg = settings->paa_segments; + const int vps = settings->ts_values_per_paa_segment; + const int ts_size = settings->timeseries_size; + const int max_card = settings->sax_alphabet_cardinality; // e.g. 256 for 8-bit SAX + if (seg <= 0 || vps <= 0 || max_card < 2) + return; + + // 1. Pool PAA values across (subsampled) series and all segments (global equi-depth). + const size_t SAMPLE_CAP = 2000000; // cap on pooled PAA values + size_t max_series = SAMPLE_CAP / (size_t)seg; + if (max_series == 0) max_series = 1; + size_t stride = 1; + if (n_series > max_series) + stride = (n_series + max_series - 1) / max_series; + + std::vector pool; + pool.reserve(((n_series + stride - 1) / stride) * (size_t)seg); + std::vector paa(seg); + for (size_t s = 0; s < n_series; s += stride) + { + paa_from_ts(database + s * (size_t)ts_size, paa.data(), seg, vps); + for (int j = 0; j < seg; j++) + pool.push_back(paa[j]); + } + if (pool.empty()) + return; + std::sort(pool.begin(), pool.end()); + + // 2. Empirical quantile function (interpolated, monotone). One shared function for + // all cardinalities preserves the iSAX nesting property (q(j/c) == q(2j/2c)). + const size_t N = pool.size(); + const float *pd = pool.data(); + auto quantile = [pd, N](double p) -> float { + if (p <= 0.0) return pd[0]; + if (p >= 1.0) return pd[N - 1]; + double x = p * (double)(N - 1); + size_t lo = (size_t)x; + if (lo + 1 >= N) return pd[N - 1]; + double frac = x - (double)lo; + return (float)(pd[lo] * (1.0 - frac) + pd[lo + 1] * frac); + }; + + // 3. Triangular table (same layout as sax_breakpoints): card c at offset (c-1)(c-2)/2. + size_t tri_size = ((size_t)max_card * (size_t)(max_card - 1)) / 2; + float *tri = (float *)malloc(sizeof(float) * tri_size); + if (tri == NULL) return; // keep Gaussian defaults on OOM + for (int c = 2; c <= max_card; c++) + { + size_t off = ((size_t)(c - 1) * (size_t)(c - 2)) / 2; + for (int j = 1; j <= c - 1; j++) + tri[off + (size_t)(j - 1)] = quantile((double)j / (double)c); + } + + // 4. Flat max-cardinality table (same layout as sax_breakpointsnew3). + int max_bp = max_card - 1; // 255 for cardinality 256 + float *flat = (float *)malloc(sizeof(float) * (size_t)max_bp); + if (flat == NULL) { free(tri); return; } + size_t off_max = ((size_t)(max_card - 1) * (size_t)(max_card - 2)) / 2; + for (int k = 0; k < max_bp; k++) + flat[k] = tri[off_max + (size_t)k]; + + // 5. Publish (freeing any previously owned tables). + if (settings->breakpoints_owned) free(settings->breakpoints_owned); + if (settings->breakpoints_max_owned) free(settings->breakpoints_max_owned); + settings->breakpoints_owned = tri; + settings->breakpoints_max_owned = flat; + settings->breakpoints = tri; + settings->breakpoints_max = flat; + } + first_buffer_layer *initialize_fbl(int initial_buffer_size, int number_of_buffers, int max_total_buffers_size, isax_index *index) { @@ -803,10 +894,10 @@ namespace daisy if (sax_value == 0 || sax_value == alphabeta_size - 1) { - return sax_breakpoints[offset + sax_value] + (sax_breakpoints[offset + sax_value] - sax_breakpoints[offset + alphabeta_size / 2 - 1]) / 2; + return daisy_active_breakpoints[offset + sax_value] + (daisy_active_breakpoints[offset + sax_value] - daisy_active_breakpoints[offset + alphabeta_size / 2 - 1]) / 2; } - return (sax_breakpoints[offset + sax_value - 1] + sax_breakpoints[offset + sax_value]) / 2; + return (daisy_active_breakpoints[offset + sax_value - 1] + daisy_active_breakpoints[offset + sax_value]) / 2; } int informed_split_decision(isax_node_split_data *split_data, @@ -872,7 +963,7 @@ namespace daisy break_point_id = (break_point_id >> ((settings->sax_bit_cardinality) - (new_bit_cardinality))) << 1; int new_cardinality = pow(2, new_bit_cardinality + 1); int offset = (new_cardinality - 1) * (new_cardinality - 2) / 2; - float b = sax_breakpoints[offset + break_point_id]; + float b = daisy_active_breakpoints[offset + break_point_id]; if (segment_to_split == -1) { @@ -1494,8 +1585,20 @@ namespace daisy int i; node->filename = (char *)malloc(sizeof(char) * index->settings->max_filename_size); - sprintf(node->filename, "%s", index->settings->root_directory); - int l = (int)strlen(index->settings->root_directory); + const char *root_dir = index->settings->root_directory; + int l; + // With no configured root directory (algorithms pass ""), contain the on-disk node + // files in a dedicated dir instead of scattering them across the CWD. + if (root_dir == NULL || root_dir[0] == '\0') + { + static const char *node_dir = "daisy_index_nodes"; + mkdir(node_dir, 0777); // idempotent + l = sprintf(node->filename, "%s/", node_dir); + } + else + { + l = sprintf(node->filename, "%s", root_dir); + } // If this has a parent then it is not a root node and as such it does have some // split data on its parent about the cardinalities. diff --git a/lib/isax/iSAXIndex.hpp b/lib/isax/iSAXIndex.hpp index 87615be..362b789 100644 --- a/lib/isax/iSAXIndex.hpp +++ b/lib/isax/iSAXIndex.hpp @@ -23,6 +23,14 @@ namespace daisy unsigned long disk_data_partial; } meminfo; + // BP_GAUSSIAN: hardcoded N(0,1) quantiles (default, for z-normalized data). + // BP_EQUIDEPTH: data-adaptive empirical quantiles computed per index at build time. + typedef enum + { + BP_GAUSSIAN = 0, + BP_EQUIDEPTH = 1 + } breakpoint_mode; + typedef struct { char new_index; @@ -65,6 +73,14 @@ namespace daisy int total_loaded_leaves; + // Gaussian: alias the global tables. Equi-depth: per-index tables (the *_owned + // pointers hold the heap allocations). breakpoints = triangular, _max = flat. + breakpoint_mode bp_mode; + const float *breakpoints; + const float *breakpoints_max; + float *breakpoints_owned; + float *breakpoints_max_owned; + } isax_index_settings; typedef struct isax_node_split_data @@ -365,7 +381,13 @@ namespace daisy int max_leaf_size, int min_leaf_size, int initial_leaf_buffer_size, int max_total_buffer_size, int initial_fbl_buffer_size, - int total_loaded_leaves, int tight_bound, int aggressive_check, int new_index, char inmemory_flag); + int total_loaded_leaves, int tight_bound, int aggressive_check, int new_index, char inmemory_flag, + breakpoint_mode bp_mode = BP_GAUSSIAN); + + // Compute per-index equi-depth breakpoints from a sample of database (row-major + // [n_series * timeseries_size]). No-op unless settings->bp_mode == BP_EQUIDEPTH. + void compute_equidepth_breakpoints(isax_index_settings *settings, + const float *database, size_t n_series); first_buffer_layer *initialize_fbl(int initial_buffer_size, int number_of_buffers, int max_total_buffers_size, isax_index *index); diff --git a/pybinds/setup.cpp b/pybinds/setup.cpp index d63507b..8dd975b 100644 --- a/pybinds/setup.cpp +++ b/pybinds/setup.cpp @@ -43,6 +43,13 @@ PYBIND11_MODULE(_core, m) .value("DTW", daisy::DistanceType::DTW) .export_values(); + ////// BREAKPOINTMODE ////// + // GAUSSIAN (default, z-normalized data) or EQUIDEPTH (data-adaptive, per index). + pybind11::enum_(m, "BreakpointMode", pybind11::arithmetic()) + .value("GAUSSIAN", daisy::BP_GAUSSIAN) + .value("EQUIDEPTH", daisy::BP_EQUIDEPTH) + .export_values(); + ////// QUERYTYPE ////// pybind11::enum_(m, "QueryType") .value("TOP_K", daisy::QueryType::TOP_K) @@ -140,6 +147,8 @@ PYBIND11_MODULE(_core, m) .def("setNumThreads", &daisy::LbBruteforce::setNumThreads, "Set the number of threads to use") .def("setPaaSegments", &daisy::LbBruteforce::setPaaSegments, "Set the number of PAA segments") .def("setSaxCardinality", &daisy::LbBruteforce::setSaxCardinality, "Set the SAX cardinality") + .def("setBreakpointMode", &daisy::LbBruteforce::setBreakpointMode, "Set breakpoint mode: 0=Gaussian (default), 1=equi-depth (data-adaptive)") + .def("getBreakpointMode", &daisy::LbBruteforce::getBreakpointMode, "Get the breakpoint mode (0=Gaussian, 1=equi-depth)") .def("setLeafSize", &daisy::LbBruteforce::setLeafSize, "Set the leaf size") .def("setMinLeafSize", &daisy::LbBruteforce::setMinLeafSize, "Set the minimum leaf size") .def("setInitialLblSize", &daisy::LbBruteforce::setInitialLblSize, "Set the initial size of the lower-bound leaf buffer (LBL)") @@ -281,6 +290,8 @@ PYBIND11_MODULE(_core, m) .def("setNumThreads", &daisy::Messi::setNumThreads, "Set the number of threads to use for both indexing and search") .def("setPaaSegments", &daisy::Messi::setPaaSegments, "Set the number of PAA segments") .def("setSaxCardinality", &daisy::Messi::setSaxCardinality, "Set the SAX cardinality") + .def("setBreakpointMode", &daisy::Messi::setBreakpointMode, "Set breakpoint mode: 0=Gaussian (default), 1=equi-depth (data-adaptive)") + .def("getBreakpointMode", &daisy::Messi::getBreakpointMode, "Get the breakpoint mode (0=Gaussian, 1=equi-depth)") .def("setLeafSize", &daisy::Messi::setLeafSize, "Set the leaf size of the index tree") .def("setMinLeafSize", &daisy::Messi::setMinLeafSize, "Set the minimum size of a leaf") .def("setInitialLblSize", &daisy::Messi::setInitialLblSize, "Set the initial LBL size") @@ -431,6 +442,8 @@ PYBIND11_MODULE(_core, m) // Setters .def("setNumThreads", &daisy::ParIS::setNumThreads, "Set the number of threads to use") + .def("setBreakpointMode", &daisy::ParIS::setBreakpointMode, "Set breakpoint mode: 0=Gaussian (default), 1=equi-depth (data-adaptive)") + .def("getBreakpointMode", &daisy::ParIS::getBreakpointMode, "Get the breakpoint mode (0=Gaussian, 1=equi-depth)") .def("setWarpingWindow", &daisy::ParIS::setWarpingWindow, "Set the warping window size for DTW (typically 10% of time series length)") // Getters @@ -573,6 +586,8 @@ PYBIND11_MODULE(_core, m) // Setters .def("setNumThreads", &daisy::Sofa::setNumThreads, "Set the number of threads to use for both indexing and search") + .def("setBreakpointMode", &daisy::Sofa::setBreakpointMode, "Set breakpoint mode: 0=Gaussian (default), 1=equi-depth (data-adaptive)") + .def("getBreakpointMode", &daisy::Sofa::getBreakpointMode, "Get the breakpoint mode (0=Gaussian, 1=equi-depth)") .def("setWordLength", &daisy::Sofa::setWordLength, "Set the SFA word length") .def("setAlphabetSize", &daisy::Sofa::setAlphabetSize, "Set the alphabet size") .def("setSearchWorkers", &daisy::Sofa::setSearchWorkers, "Set the number of search worker threads") @@ -767,6 +782,8 @@ PYBIND11_MODULE(_core, m) .def("setNumThreads", &daisy::Fresh::setNumThreads, "Set the number of threads to use for both indexing and search") .def("setPaaSegments", &daisy::Fresh::setPaaSegments, "Set the number of PAA segments") .def("setSaxCardinality", &daisy::Fresh::setSaxCardinality, "Set the SAX cardinality") + .def("setBreakpointMode", &daisy::Fresh::setBreakpointMode, "Set breakpoint mode: 0=Gaussian (default), 1=equi-depth (data-adaptive)") + .def("getBreakpointMode", &daisy::Fresh::getBreakpointMode, "Get the breakpoint mode (0=Gaussian, 1=equi-depth)") .def("setLeafSize", &daisy::Fresh::setLeafSize, "Set the leaf size of the index tree") .def("setMinLeafSize", &daisy::Fresh::setMinLeafSize, "Set the minimum size of a leaf") .def("setInitialLblSize", &daisy::Fresh::setInitialLblSize, "Set the initial LBL size")