From 837fb9ffa4e035f3bc3191ceee436658889ea840 Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Thu, 25 Jun 2026 14:11:41 +0200 Subject: [PATCH 01/19] Add DumpyOS algorithm (adapted from FADASNode reference implementation) --- demos/CMakeLists.txt | 30 ++ demos/demo_DumpyOS_L2Square.cpp | 42 +++ demos/demo_DumpyOS_L2Square.py | 36 +++ lib/algos/CMakeLists.txt | 12 + lib/algos/DumpyOS.cpp | 533 ++++++++++++++++++++++++++++++++ lib/algos/DumpyOS.hpp | 52 ++++ lib/daisy.hpp | 1 + pybinds/setup.cpp | 40 +++ 8 files changed, 746 insertions(+) create mode 100644 demos/demo_DumpyOS_L2Square.cpp create mode 100644 demos/demo_DumpyOS_L2Square.py create mode 100644 lib/algos/DumpyOS.cpp create mode 100644 lib/algos/DumpyOS.hpp diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt index dc71dae..6e330d7 100644 --- a/demos/CMakeLists.txt +++ b/demos/CMakeLists.txt @@ -436,6 +436,36 @@ if(BUILD_DEMO) message(STATUS "Include directories added for demo_Hercules_L2Square.") endif() + # ////// DUMPYOS L2Square ////// + if(DEBUG_MSG) + message(STATUS "---") + message(STATUS "## Demo: DumpyOS L2Square") + message(STATUS "Attempting to add executable: demo_DumpyOS_L2Square") + endif() + add_executable(demo_DumpyOS_L2Square demo_DumpyOS_L2Square.cpp) + if(DEBUG_MSG) + message(STATUS "Executable demo_DumpyOS_L2Square added.") + endif() + + if(DEBUG_MSG) + message(STATUS "Linking libraries for demo_DumpyOS_L2Square...") + endif() + target_link_libraries(demo_DumpyOS_L2Square PRIVATE dino_lib commons_lib) + if(DEBUG_MSG) + message(STATUS "Libraries linked for demo_DumpyOS_L2Square.") + endif() + + if(DEBUG_MSG) + message(STATUS "Adding include directories for demo_DumpyOS_L2Square...") + endif() + target_include_directories(demo_DumpyOS_L2Square PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + if(DEBUG_MSG) + message(STATUS "Include directories added for demo_DumpyOS_L2Square.") + endif() + else() if(DEBUG_MSG) message(STATUS "BUILD_DEMO is FALSE. Demo executables will NOT be built.") diff --git a/demos/demo_DumpyOS_L2Square.cpp b/demos/demo_DumpyOS_L2Square.cpp new file mode 100644 index 0000000..2c4d6a4 --- /dev/null +++ b/demos/demo_DumpyOS_L2Square.cpp @@ -0,0 +1,42 @@ +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" +#include + +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::DumpyOS dumpyos_search(daisy::DistanceType::L2_SQUARED); + dumpyos_search.setNumThreads(4); + + dumpyos_search.buildIndex(database, n_database, dim); + + daisy::idx_t *I = new daisy::idx_t[n_query * k]; + float *D = new float[n_query * k]; + dumpyos_search.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_DumpyOS_L2Square.py b/demos/demo_DumpyOS_L2Square.py new file mode 100644 index 0000000..842aafd --- /dev/null +++ b/demos/demo_DumpyOS_L2Square.py @@ -0,0 +1,36 @@ +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, DumpyOS + +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 = DumpyOS(DistanceType.L2_SQUARED) + + index.setNumThreads(4) + index.buildIndex(db) + + 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/lib/algos/CMakeLists.txt b/lib/algos/CMakeLists.txt index 56ef9c8..33bc218 100644 --- a/lib/algos/CMakeLists.txt +++ b/lib/algos/CMakeLists.txt @@ -152,6 +152,18 @@ if(DEBUG_MSG) message(STATUS "Hercules.cpp added.") endif() +# ////// DUMPYOS ////// +if(DEBUG_MSG) + message(STATUS "Adding DumpyOS.cpp to dino_lib sources.") +endif() +target_sources(dino_lib + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/DumpyOS.cpp +) +if(DEBUG_MSG) + message(STATUS "DumpyOS.cpp added.") +endif() + # ////// PARIS ////// if(DEBUG_MSG) message(STATUS "Adding ParIS.cpp to dino_lib sources.") diff --git a/lib/algos/DumpyOS.cpp b/lib/algos/DumpyOS.cpp new file mode 100644 index 0000000..2678a95 --- /dev/null +++ b/lib/algos/DumpyOS.cpp @@ -0,0 +1,533 @@ +// Adapted from FADASNode.cpp in DumpyOS/src/IndexConstruction/FADASNode.cpp +// (Wang et al., "DumpyOS: A Dumpy Index for Scalable Data Series Similarity +// Search", VLDB Journal 2024) +// Reference: SaxUtil.cpp (getMidLineFromSaxSymbolbc8, extendSax), +// MathUtil.cpp (deviation, generateMaskSettingKbits) + +#include "DumpyOS.hpp" +#include "../isax/SAX.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace daisy { + +// ---- bp8 : normal-distribution breakpoints for 8-bit SAX ---- +// Directly from SaxUtil.cpp in the DumpyOS reference implementation. +// bp8[i] is the upper boundary of the i-th SAX symbol region. +// bp8[255] = +inf sentinel; midpoints are extrapolated for symbols 0 and 255. +static const double bp8[256] = { + -2.660067468617458,-2.4175590162365035,-2.2662268092096522, + -2.1538746940614573,-2.063527898316245,-1.9874278859298962, + -1.921350774293703,-1.8627318674216515,-1.8098922384806087, + -1.7616704103630665,-1.7172281175057411,-1.6759397227734438, + -1.637325382768064,-1.601008664886076,-1.5666885860684134, + -1.534120544352546,-1.5031029431292737,-1.4734675779471014, + -1.4450725798180746,-1.4177971379962677,-1.3915374879959008, + -1.3662038163720984,-1.341717841080254,-1.318010897303537, + -1.2950224067058147,-1.2726986411905359,-1.2509917154625454, + -1.229858759216589,-1.209261231709155,-1.189164350199337, + -1.169536610207143,-1.1503493803760083,-1.131576558386188, + -1.113194277160929,-1.0951806527613883,-1.0775155670402805, + -1.0601804794353549,-1.0431582633184537,-1.0264330631379108, + -1.0099901692495823,-0.993815907860883,-0.9778975439405418, + -0.9622231952954206,-0.946781756301046,-0.9315628300071148, + -0.9165566675331128,-0.9017541138301002,-0.8871465590188762, + -0.8727258946270402,-0.858484474141832,-0.8444150773752572, + -0.8305108782053992,-0.8167654153150912,-0.8031725655979178, + -0.7897265199432658,-0.7764217611479275,-0.7632530437325706, + -0.7502153754679404,-0.7373040004386545,-0.7245143834923653, + -0.711842195939419,-0.69928330238322,-0.6868337485747303, + -0.6744897501960819,-0.6622476824884141,-0.6501040706479954, + -0.6380555809225171,-0.6260990123464211,-0.6142312890602454, + -0.6024494531644237,-0.5907506580628189,-0.5791321622555561, + -0.5675913235445692,-0.5561255936186916,-0.5447325129881759, + -0.5334097062412806,-0.5221548775980015,-0.5109658067382474, + -0.4998403448837353,-0.4887764111146696,-0.4777719889038861, + -0.46682512285258965,-0.4559339156131388,-0.44509652498551644, + -0.4343111611752096,-0.42357608420119974,-0.41288960144365433, + -0.40225006532172536,-0.3916558710925915,-0.3811054547635565, + -0.3705972911096293,-0.36012989178956945,-0.3497018035538953, + -0.3393116065388173,-0.3289579126404911,-0.31863936396437526, + -0.30835463134483726,-0.2981024129304869,-0.2878814328310118, + -0.27769043982157676,-0.2675282061010972,-0.25739352610093835, + -0.24728521534080486,-0.2372021093287877,-0.22714306250271535, + -0.2171069472101298,-0.2070926527243603,-0.19709908429431236, + -0.18712516222572084,-0.17716982099173986,-0.16723200837085014, + -0.15731068461017073,-0.14740482161235488,-0.13751340214433597, + -0.12763541906627035,-0.11776987457909531,-0.10791577948918657, + -0.0980721524886611,-0.08823801944992447,-0.07841241273311222, + -0.06859437050511813,-0.05878293606894307,-0.04897715720213194, + -0.03917608550309764,-0.02937877574415705,-0.019584285230126924, + -0.009791673161345348, + 0.0, + 0.009791673161345348, 0.019584285230126924, 0.02937877574415705, + 0.03917608550309764, 0.04897715720213194, 0.05878293606894307, + 0.06859437050511813, 0.07841241273311222, 0.08823801944992447, + 0.0980721524886611, 0.10791577948918657, 0.11776987457909531, + 0.12763541906627035, 0.13751340214433597, 0.14740482161235488, + 0.15731068461017073, 0.16723200837085014, 0.17716982099173986, + 0.18712516222572084, 0.19709908429431236, 0.2070926527243603, + 0.2171069472101298, 0.22714306250271535, 0.2372021093287877, + 0.24728521534080486, 0.25739352610093835, 0.2675282061010972, + 0.27769043982157676, 0.2878814328310118, 0.2981024129304869, + 0.30835463134483726, 0.31863936396437526, 0.3289579126404911, + 0.3393116065388173, 0.3497018035538953, 0.36012989178956945, + 0.3705972911096293, 0.3811054547635565, 0.3916558710925915, + 0.40225006532172536, 0.41288960144365433, 0.42357608420119974, + 0.4343111611752096, 0.44509652498551644, 0.4559339156131388, + 0.46682512285258965, 0.4777719889038861, 0.4887764111146696, + 0.4998403448837353, 0.5109658067382474, 0.5221548775980015, + 0.5334097062412806, 0.5447325129881759, 0.5561255936186916, + 0.5675913235445692, 0.5791321622555561, 0.5907506580628189, + 0.6024494531644237, 0.6142312890602454, 0.6260990123464211, + 0.6380555809225171, 0.6501040706479954, 0.6622476824884141, + 0.6744897501960819, 0.6868337485747303, 0.69928330238322, + 0.711842195939419, 0.7245143834923653, 0.7373040004386545, + 0.7502153754679404, 0.7632530437325706, 0.7764217611479275, + 0.7897265199432658, 0.8031725655979178, 0.8167654153150912, + 0.8305108782053992, 0.8444150773752572, 0.858484474141832, + 0.8727258946270402, 0.8871465590188762, 0.9017541138301002, + 0.9165566675331128, 0.9315628300071148, 0.946781756301046, + 0.9622231952954206, 0.9778975439405418, 0.993815907860883, + 1.0099901692495823, 1.0264330631379108, 1.0431582633184537, + 1.0601804794353549, 1.0775155670402805, 1.0951806527613883, + 1.113194277160929, 1.131576558386188, 1.1503493803760083, + 1.169536610207143, 1.189164350199337, 1.209261231709155, + 1.229858759216589, 1.2509917154625454, 1.2726986411905359, + 1.2950224067058147, 1.318010897303537, 1.341717841080254, + 1.3662038163720984, 1.3915374879959008, 1.4177971379962677, + 1.4450725798180746, 1.4734675779471014, 1.5031029431292737, + 1.534120544352546, 1.5666885860684134, 1.601008664886076, + 1.637325382768064, 1.6759397227734438, 1.7172281175057411, + 1.7616704103630665, 1.8098922384806087, 1.8627318674216515, + 1.921350774293703, 1.9874278859298962, 2.063527898316245, + 2.1538746940614573, 2.2662268092096522, 2.4175590162365035, + 2.660067468617458, + // sentinel for symbol 255 upper bound (numeric_limits::max()) + 3.4028234663852886e+38 +}; + +// ---- internal helpers (adapted from SaxUtil and MathUtil) ---- + +// Equivalent to SaxUtil::getMidLineFromSaxSymbolbc8. +static double get_midpoint(unsigned sym) { + if (sym == 0) return bp8[0] - (bp8[1] - bp8[0]); + if (sym == 255) return bp8[254] + (bp8[254] - bp8[253]); + return (bp8[sym - 1] + bp8[sym]) * 0.5; +} + +// Equivalent to SaxUtil::extendSax(sax, bits_cardinality) — all segments. +// Produces a vertexNum (1<= 0) ? ((sax[s] >> shift) & 1) : 0; + res = (res << 1) | bit; + } + return res; +} + +// Equivalent to SaxUtil::extendSax(sax, bits_cardinality, chosenSegments). +// Produces a child ID in [0, 2^|chosen|). +static int extend_sax_chosen(const sax_type* sax, const int* levels, + const std::vector& chosen, int max_bits) { + int res = 0; + for (int s : chosen) { + int shift = max_bits - levels[s] - 1; + int bit = (shift >= 0) ? ((sax[s] >> shift) & 1) : 0; + res = (res << 1) | bit; + } + return res; +} + +// Equivalent to MathUtil::generateMaskSettingKbits(plan, k, len). +// plan must be sorted ascending. Sets bit (w-1-plan[i]) for each i in plan. +static int generate_mask(const int* plan, int lambda, int w) { + int res = 0, cur = 0; + for (int i = 0; i < w; ++i) { + res <<= 1; + if (cur < lambda && plan[cur] == i) { res |= 1; ++cur; } + } + return res; +} + +// Population standard deviation — equivalent to MathUtil::deviation(double*, int). +static double pop_stdev(const double* vals, int n) { + double mean = 0.0; + for (int i = 0; i < n; ++i) mean += vals[i]; + mean /= n; + double var = 0.0; + for (int i = 0; i < n; ++i) { double d = vals[i] - mean; var += d * d; } + return std::sqrt(var / n); +} + +// Equivalent to FADASNode::compute_score. +// node_n = node->n (total series in node) +// alpha = config_.alpha +static double compute_score_(const std::vector& node_sizes, + const int* plan, int lambda, + const std::vector& data_seg_stdev, + int node_n, int leaf_size, double alpha) { + if (node_n < 2 * leaf_size) { + if ((int)node_sizes.size() >= 2 && + (node_sizes[0] > leaf_size || node_sizes[1] > leaf_size)) + return (double)std::min(node_sizes[0], node_sizes[1]) / leaf_size; + return data_seg_stdev[plan[0]] * 100.0; + } + int over = 0; + for (int s : node_sizes) if (s > leaf_size) ++over; + double w_frac = (double)over / (int)node_sizes.size(); + + double sum_seg = 0.0; + for (int i = 0; i < lambda; ++i) sum_seg += data_seg_stdev[plan[i]]; + sum_seg = std::exp(1.0 + std::sqrt(sum_seg / lambda)); + + int nc = (int)node_sizes.size(); + std::vector tmp(nc); + for (int i = 0; i < nc; ++i) tmp[i] = (double)node_sizes[i] / leaf_size; + double sigma_f = pop_stdev(tmp.data(), nc); + double balance = std::exp(-(1.0 + w_frac) * sigma_f); + return sum_seg + alpha * balance; +} + +// Equivalent to FADASNode::visitPlanFromBaseTable. +// cur_lambda = parent lambda minus 1 (we're evaluating plans of size cur_lambda). +// plan = parent plan of size cur_lambda + 1 (sorted ascending). +// base_tbl = parent child fill counts of size 1 << (cur_lambda + 1). +static void visit_plan_from_base_table( + std::unordered_set& visited, + int cur_lambda, const int* plan, const std::vector& base_tbl, + double* max_score, std::vector& best_plan, + int lambda_min, int mask_code, + const std::vector& data_seg_stdev, + int node_n, int leaf_size, double alpha, int w) +{ + // base_mask has (cur_lambda+1) bits set (covers all 1<<(cur_lambda+1) entries) + int base_mask = 1; + for (int i = 0; i < cur_lambda; ++i) base_mask = (base_mask << 1) | 1; + + for (int i = 0; i <= cur_lambda; ++i) { + int reset_pos = plan[i]; + int cur_whole_mask = mask_code - (1 << (w - 1 - reset_pos)); + if (visited.count(cur_whole_mask)) continue; + visited.insert(cur_whole_mask); + + // Sub-plan: drop element i from parent plan + int* new_plan = new int[cur_lambda]; + for (int j = 0, k = 0; j <= cur_lambda; ++j) + if (j != i) new_plan[k++] = plan[j]; + + // Aggregate base_tbl into cur_lambda child fill counts by masking out bit i + int cur_base_mask = base_mask - (1 << (cur_lambda - i)); + std::map nsmap; + for (int j = 0; j < (int)base_tbl.size(); ++j) + nsmap[cur_base_mask & j] += base_tbl[j]; + std::vector new_tbl; + new_tbl.reserve(1 << cur_lambda); + for (auto& kv : nsmap) new_tbl.push_back(kv.second); + + double score = compute_score_(new_tbl, new_plan, cur_lambda, + data_seg_stdev, node_n, leaf_size, alpha); + if (score > *max_score) { + *max_score = score; + best_plan.assign(new_plan, new_plan + cur_lambda); + } + if (cur_lambda > lambda_min) + visit_plan_from_base_table(visited, cur_lambda - 1, new_plan, new_tbl, + max_score, best_plan, lambda_min, + cur_whole_mask, data_seg_stdev, + node_n, leaf_size, alpha, w); + delete[] new_plan; + } +} + +// Equivalent to FADASNode::determineFanout. +static void determine_fanout(int n, int leaf_size, double f_low, double f_high, int w, + int* lambda_min, int* lambda_max) { + if (n < 2 * leaf_size) { *lambda_min = 1; *lambda_max = 1; return; } + *lambda_min = -1; + *lambda_max = w; + int vertex_num = 1 << w; + double _min = (double)n / (leaf_size * f_high); + double _max = (double)n / (leaf_size * f_low); + if ((double)vertex_num < _min) { *lambda_min = w; *lambda_max = w; return; } + for (int i = 1; i <= w; ++i) { + if (*lambda_min == -1) { + if ((double)(1 << i) >= _min) *lambda_min = i; + } else { + if ((double)(1 << i) == _max) { *lambda_max = i; break; } + else if ((double)(1 << i) > _max) { *lambda_max = std::max(i - 1, *lambda_min); break; } + } + } + if (*lambda_min == -1) { *lambda_min = w; *lambda_max = w; } +} + +// Equivalent to FADASNode::determineSegments. +void DumpyOS::determineSegments_(DumpyOSNode* node) { + int w = config_.paa_segments; + int max_bits = config_.sax_bit_cardinality; + int ls = config_.leaf_size; + double alpha = config_.alpha; + + int lambda_min, lambda_max; + determine_fanout(node->n, ls, config_.fill_lower, config_.fill_upper, + w, &lambda_min, &lambda_max); + + // Edge case: need all segments (vertexNum < _min) + if (lambda_min == w && lambda_max == w) { + for (int i = 0; i < w; ++i) node->chosen_segs.push_back(i); + return; + } + + int vertex_num = 1 << w; + + // Build unit_size table: aggregate series counts into 2^w first-level buckets + std::vector unit_size(vertex_num, 0); + // Per-segment SAX symbol frequency (for data_seg_stdev) + std::vector> seg_sym_cnt(w); + + for (idx_t si : node->entries) { + const sax_type* sax = sax_table_ + (size_t)si * w; + for (int i = 0; i < w; ++i) + seg_sym_cnt[i][(int)(unsigned char)sax[i]]++; + int head = extend_sax_all(sax, node->levels.data(), w, max_bits); + unit_size[head]++; + } + + // Compute data_seg_mean and data_seg_stdev (variance of midpoint values) + // Equivalent to the data_seg_mean/stdev loop in FADASNode::determineSegments. + std::vector data_seg_mean(w, 0.0), data_seg_stdev(w, 0.0); + for (int i = 0; i < w; ++i) { + for (auto& kv : seg_sym_cnt[i]) + data_seg_mean[i] += get_midpoint((unsigned)kv.first) * kv.second; + data_seg_mean[i] /= node->n; + for (auto& kv : seg_sym_cnt[i]) { + double d = get_midpoint((unsigned)kv.first) - data_seg_mean[i]; + data_seg_stdev[i] += kv.second * d * d; + } + data_seg_stdev[i] /= node->n; + } + + double max_score = 0.0; + std::vector best_plan; + std::unordered_set visited; + + // Enumerate all C(w, lambda_max) plans; for each, also search sub-plans + // via visitPlanFromBaseTable (lambda_min..lambda_max-1). + std::vector idx(lambda_max); + std::iota(idx.begin(), idx.end(), 0); + + while (true) { + const int* plan = idx.data(); + + // Compute child fill counts using unit_size table + bit masking + int mask_code = generate_mask(plan, lambda_max, w); + std::map node_size_map; + for (int j = 0; j < vertex_num; ++j) + node_size_map[mask_code & j] += unit_size[j]; + std::vector plan_node_sizes; + plan_node_sizes.reserve(1 << lambda_max); + for (auto& kv : node_size_map) plan_node_sizes.push_back(kv.second); + + double score = compute_score_(plan_node_sizes, plan, lambda_max, + data_seg_stdev, node->n, ls, alpha); + if (score > max_score) { + max_score = score; + best_plan.assign(plan, plan + lambda_max); + } + + if (lambda_min <= lambda_max - 1) + visit_plan_from_base_table(visited, lambda_max - 1, plan, plan_node_sizes, + &max_score, best_plan, lambda_min, mask_code, + data_seg_stdev, node->n, ls, alpha, w); + + // Advance to next C(w, lambda_max) combination (Gosper-style) + int i = lambda_max - 1; + while (i >= 0 && idx[i] == w - lambda_max + i) --i; + if (i < 0) break; + ++idx[i]; + for (int j = i + 1; j < lambda_max; ++j) idx[j] = idx[j - 1] + 1; + } + + node->chosen_segs = best_plan; +} + +// ---- splitNode_ ---- + +void DumpyOS::splitNode_(DumpyOSNode* node) { + int w = config_.paa_segments; + int max_bits = config_.sax_bit_cardinality; + + determineSegments_(node); + if (node->chosen_segs.empty()) return; // exhausted bit budget + + int lam = (int)node->chosen_segs.size(); + int num_ch = 1 << lam; + node->children.assign(num_ch, nullptr); + + // Create all children with inherited levels (+ 1 for each chosen segment) + for (int sid = 0; sid < num_ch; ++sid) { + DumpyOSNode* child = new DumpyOSNode(); + child->levels = node->levels; + for (int cs : node->chosen_segs) child->levels[cs]++; + node->children[sid] = child; + } + + // Redistribute series — equivalent to FADASNode's child routing + for (idx_t si : node->entries) { + const sax_type* sax = sax_table_ + (size_t)si * w; + int sid = extend_sax_chosen(sax, node->levels.data(), node->chosen_segs, max_bits); + node->children[sid]->entries.push_back(si); + node->children[sid]->n++; + } + node->entries.clear(); + node->entries.shrink_to_fit(); + + for (DumpyOSNode* child : node->children) + if (child && child->n > config_.leaf_size) + splitNode_(child); +} + +// ---- buildIndex ---- + +void DumpyOS::buildIndex(DataSource* data_source) { + dim = data_source->getDim(); + n_database = data_source->getTotalRecords(); + + const float* raw = data_source->rawPointer(); + if (raw) { + database = const_cast(raw); + owns_database_ = false; + } else { + database = new float[(size_t)n_database * dim]; + owns_database_ = true; + float* ptr = database; + while (data_source->nextRecord(ptr)) ptr += dim; + } + + int w = config_.paa_segments; + int max_bits = config_.sax_bit_cardinality; + int cardinality = 1 << max_bits; + int pts_per_seg = (int)dim / w; + + // Compute and cache the SAX word for every database series + sax_table_ = new sax_type[(size_t)n_database * w]; + std::vector paa(w); + + for (idx_t i = 0; i < n_database; ++i) { + paa_from_ts(database + (size_t)i * dim, paa.data(), w, pts_per_seg); + sax_from_paa(paa.data(), sax_table_ + (size_t)i * w, w, cardinality, max_bits); + } + + // Root: all series, all segment levels at 0 + root_ = new DumpyOSNode(); + root_->levels.assign(w, 0); + root_->n = (int)n_database; + root_->entries.resize(n_database); + std::iota(root_->entries.begin(), root_->entries.end(), (idx_t)0); + + if (root_->n > config_.leaf_size) + splitNode_(root_); +} + +// ---- searchIndex ---- + +static float l2sq_early(const float* a, const float* b, int d, float bound) { + float s = 0.0f; + for (int i = 0; i < d; ++i) { + float v = a[i] - b[i]; + s += v * v; + if (s >= bound) return s; + } + return s; +} + +void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, + idx_t* I, float* D) { + if (!validateSearchParams(k, n_query)) return; + + int w = config_.paa_segments; + int max_bits = config_.sax_bit_cardinality; + int cardinality = 1 << max_bits; + int pts_per_seg = (int)dim / w; + + std::vector q_sax(w); + std::vector q_paa(w); + + for (idx_t qi = 0; qi < n_query; ++qi) { + const float* q = query + (size_t)qi * dim; + + paa_from_ts(q, q_paa.data(), w, pts_per_seg); + sax_from_paa(q_paa.data(), q_sax.data(), w, cardinality, max_bits); + + // Route to leaf — equivalent to FADASNode::route() + DumpyOSNode* node = root_; + while (!node->chosen_segs.empty()) { + int sid = extend_sax_chosen(q_sax.data(), node->levels.data(), + node->chosen_segs, max_bits); + if (sid < 0 || sid >= (int)node->children.size() || + node->children[sid] == nullptr) break; + node = node->children[sid]; + } + + // kNN via max-heap + using Pair = std::pair; + std::vector heap; + heap.reserve((size_t)k + 1); + float bsf = FLT_MAX; + + for (idx_t si : node->entries) { + float dist = l2sq_early(q, database + (size_t)si * dim, (int)dim, bsf); + if ((idx_t)heap.size() < k || dist < bsf) { + heap.push_back({dist, si}); + std::push_heap(heap.begin(), heap.end()); + if ((idx_t)heap.size() > k) { + std::pop_heap(heap.begin(), heap.end()); + heap.pop_back(); + } + if ((idx_t)heap.size() == k) bsf = heap.front().first; + } + } + + std::sort_heap(heap.begin(), heap.end()); + for (idx_t j = 0; j < k; ++j) { + if (j < (idx_t)heap.size()) { + I[qi * k + j] = heap[j].second; + D[qi * k + j] = heap[j].first; + } else { + I[qi * k + j] = 0; + D[qi * k + j] = FLT_MAX; + } + } + } +} + +// ---- destructor ---- + +void DumpyOS::destroyTree_(DumpyOSNode* node) { + if (!node) return; + for (DumpyOSNode* child : node->children) destroyTree_(child); + delete node; +} + +DumpyOS::DumpyOS(DistanceType distance_type) + : SimilaritySearchAlgorithm(distance_type) {} + +DumpyOS::DumpyOS(DistanceType distance_type, const DumpyOSConfig& config) + : SimilaritySearchAlgorithm(distance_type), config_(config) {} + +DumpyOS::~DumpyOS() { + destroyTree_(root_); + delete[] sax_table_; + if (owns_database_) delete[] database; +} + +} // namespace daisy diff --git a/lib/algos/DumpyOS.hpp b/lib/algos/DumpyOS.hpp new file mode 100644 index 0000000..8c003e8 --- /dev/null +++ b/lib/algos/DumpyOS.hpp @@ -0,0 +1,52 @@ +#pragma once + +#include "SimilaritySearchAlgorithm.hpp" +#include "../isax/iSAXTypes.hpp" + +#include + +namespace daisy { + +struct DumpyOSConfig { + int leaf_size = 10000; + int paa_segments = 16; + int sax_bit_cardinality = 8; + float alpha = 0.2f; + float fill_lower = 0.5f; // f_low in the paper + float fill_upper = 3.0f; // f_high in the paper +}; + +// Adapted from FADASNode in DumpyOS/include/DataStructures/FADASNode.h +struct DumpyOSNode { + std::vector levels; // bits_cardinality[] per segment + std::vector chosen_segs; // chosenSegments; empty ↔ leaf + std::vector children; // 2^|chosen_segs| entries (may be nullptr) + std::vector entries; // series indices (leaf only) + int n = 0; +}; + +class DumpyOS : public SimilaritySearchAlgorithm { +public: + DumpyOS(DistanceType distance_type); + DumpyOS(DistanceType distance_type, const DumpyOSConfig& config); + + using SimilaritySearchAlgorithm::buildIndex; + + void buildIndex(DataSource* data_source) override; + void searchIndex(const float* query, idx_t n_query, idx_t k, + idx_t* I, float* D) override; + + ~DumpyOS() override; + +private: + DumpyOSConfig config_; + DumpyOSNode* root_ = nullptr; + sax_type* sax_table_ = nullptr; // [n_database * paa_segments] + bool owns_database_ = false; + + void determineSegments_(DumpyOSNode* node); + void splitNode_(DumpyOSNode* node); + void destroyTree_(DumpyOSNode* node); +}; + +} // namespace daisy diff --git a/lib/daisy.hpp b/lib/daisy.hpp index 5c4411f..d34ac2c 100644 --- a/lib/daisy.hpp +++ b/lib/daisy.hpp @@ -9,5 +9,6 @@ #include "algos/hodyssey/Odyssey.hpp" #include "algos/Sofa.hpp" #include "algos/Hercules.hpp" +#include "algos/DumpyOS.hpp" #endif diff --git a/pybinds/setup.cpp b/pybinds/setup.cpp index a39c8d4..634b671 100644 --- a/pybinds/setup.cpp +++ b/pybinds/setup.cpp @@ -24,6 +24,7 @@ #endif #include "../lib/algos/DataSource.hpp" #include "../lib/algos/Hercules.hpp" +#include "../lib/algos/DumpyOS.hpp" #ifdef SOFA_FFTW_ENABLED #if SOFA_FFTW_ENABLED != 0 #include "../lib/algos/Sofa.hpp" @@ -597,4 +598,43 @@ PYBIND11_MODULE(_core, m) pybind11::array_t({n_query, k}, indices.data()), pybind11::array_t({n_query, k}, distances.data()) ); }, "Search the Hercules index and return (indices, distances)"); + + pybind11::class_(m, "DumpyOSConfig", "Configuration for the DumpyOS similarity search index") + .def(pybind11::init<>()) + .def_readwrite("leaf_size", &daisy::DumpyOSConfig::leaf_size) + .def_readwrite("paa_segments", &daisy::DumpyOSConfig::paa_segments) + .def_readwrite("sax_bit_cardinality", &daisy::DumpyOSConfig::sax_bit_cardinality) + .def_readwrite("alpha", &daisy::DumpyOSConfig::alpha) + .def_readwrite("fill_lower", &daisy::DumpyOSConfig::fill_lower) + .def_readwrite("fill_upper", &daisy::DumpyOSConfig::fill_upper); + + pybind11::class_(m, "DumpyOS", "DumpyOS iSAX-based multi-ary adaptive time series similarity index") + .def(pybind11::init(), "Create a new DumpyOS with the given distance metric") + .def(pybind11::init(), "Create a new DumpyOS with the given distance metric and configuration") + .def("setNumThreads", &daisy::DumpyOS::setNumThreads, "Set the number of threads") + .def("buildIndex", [](daisy::DumpyOS &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 DumpyOS index from a 2D float32 NumPy array") + .def("searchIndex", [](daisy::DumpyOS &self, pybind11::array_t query, daisy::idx_t k) + { + pybind11::buffer_info query_buf = query.request(); + if (query_buf.ndim != 2) + throw std::runtime_error("Query array must be 2D"); + if (k <= 0) + throw std::runtime_error("k must be positive"); + const daisy::idx_t n_query = query_buf.shape[0]; + std::vector indices(n_query * k); + std::vector distances(n_query * k); + self.searchIndex(static_cast(query_buf.ptr), n_query, k, + indices.data(), distances.data()); + return pybind11::make_tuple( + pybind11::array_t({n_query, k}, indices.data()), + pybind11::array_t({n_query, k}, distances.data()) + ); }, "Search the DumpyOS index and return (indices, distances)"); } From 143fc9904f586b97a8caf6041ce50a4829254cbe Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Thu, 25 Jun 2026 16:30:17 +0200 Subject: [PATCH 02/19] dumpy tests --- tests/CMakeLists.txt | 27 ++++++++++++++++++++++ tests/test_DumpyOS_L2Square.cpp | 41 +++++++++++++++++++++++++++++++++ tests/test_utils.hpp | 11 +++++++++ 3 files changed, 79 insertions(+) create mode 100644 tests/test_DumpyOS_L2Square.cpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c44061c..10eb86a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -624,4 +624,31 @@ target_include_directories(test_Hercules_L2Square gtest_discover_tests( test_Hercules_L2Square WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} +) + +# ////// DUMPYOS L2Square ////// +add_executable( + test_DumpyOS_L2Square + test_DumpyOS_L2Square.cpp + test_utils.cpp +) + +target_link_libraries( + test_DumpyOS_L2Square + PRIVATE + GTest::gtest_main + dino_lib + commons_lib + stdc++fs +) + +target_include_directories(test_DumpyOS_L2Square + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons +) + +gtest_discover_tests( + test_DumpyOS_L2Square + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} ) \ No newline at end of file diff --git a/tests/test_DumpyOS_L2Square.cpp b/tests/test_DumpyOS_L2Square.cpp new file mode 100644 index 0000000..3fb0c01 --- /dev/null +++ b/tests/test_DumpyOS_L2Square.cpp @@ -0,0 +1,41 @@ +#include "test_utils.hpp" +#include "../commons/test_bm_utils.hpp" +#include "../commons/paramSetup.hpp" + +std::string prefix = "bruteForce"; + +TEST_P(DumpyOSParameterizedTest, AllConfigurations) +{ + const SSTestConfig &config = GetParam(); + + daisy::DumpyOS search(daisy::DistanceType::L2_SQUARED); + + std::string gt_I_path = config.gt_I_prefix + std::to_string(config.k_value) + ".txt"; + std::string gt_D_path = config.gt_D_prefix + std::to_string(config.k_value) + ".txt"; + + runSST( + &search, + prefix, + gt_I_path, + gt_D_path, + config.dataset_path, + config.query_path, + config.thread_count); +} + +INSTANTIATE_TEST_SUITE_P( + DumpyOSTests, + DumpyOSParameterizedTest, + ::testing::ValuesIn(test_configs), + [](const ::testing::TestParamInfo &info) + { + return info.param.name + "_k" + std::to_string(info.param.k_value) + + "_thread" + std::to_string(info.param.thread_count) + + "_idx" + std::to_string(info.index); + }); + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/tests/test_utils.hpp b/tests/test_utils.hpp index 0a2c563..3b52b2b 100644 --- a/tests/test_utils.hpp +++ b/tests/test_utils.hpp @@ -15,6 +15,7 @@ #include "../lib/algos/Sing.hpp" #include "../lib/algos/Sofa.hpp" #include "../lib/algos/Hercules.hpp" +#include "../lib/algos/DumpyOS.hpp" class SimilaritySearchTest : public ::testing::Test { @@ -142,4 +143,14 @@ class HerculesParameterizedTest : public SimilaritySearchTest, static void TearDownTestSuite() {} }; +class DumpyOSParameterizedTest : public SimilaritySearchTest, + public ::testing::WithParamInterface +{ +protected: + using SimilaritySearchTest::runSST; + + static void SetUpTestSuite() {} + static void TearDownTestSuite() {} +}; + #endif From 413e09df597269084bc7606530d0b741c7a9c315 Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Fri, 26 Jun 2026 10:44:07 +0200 Subject: [PATCH 03/19] Fix exact search, implement pybinds and tests + updated readme --- benchmark/CMakeLists.txt | 21 +++++ lib/algos/DumpyOS.cpp | 165 +++++++++++++++++++++++++++++++-------- lib/algos/DumpyOS.hpp | 1 + 3 files changed, 156 insertions(+), 31 deletions(-) diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt index 9da0d89..a2eb474 100644 --- a/benchmark/CMakeLists.txt +++ b/benchmark/CMakeLists.txt @@ -261,6 +261,27 @@ target_include_directories(bm_Sofa_L2Square ) endif() +# ////// DUMPYOS ////// +add_executable(bm_DumpyOS_L2Square + bm_DumpyOS_L2Square.cpp + bm_utils.cpp + ../commons/paramSetup.cpp + ../commons/test_bm_utils.cpp + ../commons/dataloaders.cpp +) +target_link_libraries(bm_DumpyOS_L2Square + PRIVATE + benchmark::benchmark + benchmark::benchmark_main + dino_lib + GTest::gtest +) +target_include_directories(bm_DumpyOS_L2Square + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons +) + # ////// HERCULES ////// add_executable(bm_Hercules_L2Square bm_Hercules_L2Square.cpp diff --git a/lib/algos/DumpyOS.cpp b/lib/algos/DumpyOS.cpp index 2678a95..5e740e9 100644 --- a/lib/algos/DumpyOS.cpp +++ b/lib/algos/DumpyOS.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include @@ -373,11 +374,19 @@ void DumpyOS::splitNode_(DumpyOSNode* node) { int num_ch = 1 << lam; node->children.assign(num_ch, nullptr); - // Create all children with inherited levels (+ 1 for each chosen segment) + // Create all children with inherited levels/sax_word (+ 1 bit for each chosen segment) for (int sid = 0; sid < num_ch; ++sid) { DumpyOSNode* child = new DumpyOSNode(); - child->levels = node->levels; - for (int cs : node->chosen_segs) child->levels[cs]++; + child->levels = node->levels; + child->sax_word = node->sax_word; + int cid = sid; + // Extract bits LSB-first: last chosen_seg gets bit 0 of sid + for (int j = lam - 1; j >= 0; --j) { + int cs = node->chosen_segs[j]; + child->sax_word[cs] = (child->sax_word[cs] << 1) | (cid & 1); + child->levels[cs]++; + cid >>= 1; + } node->children[sid] = child; } @@ -427,9 +436,10 @@ void DumpyOS::buildIndex(DataSource* data_source) { sax_from_paa(paa.data(), sax_table_ + (size_t)i * w, w, cardinality, max_bits); } - // Root: all series, all segment levels at 0 + // Root: all series, all segment levels and sax_word at 0 root_ = new DumpyOSNode(); root_->levels.assign(w, 0); + root_->sax_word.assign(w, 0); root_->n = (int)n_database; root_->entries.resize(n_database); std::iota(root_->entries.begin(), root_->entries.end(), (idx_t)0); @@ -440,6 +450,53 @@ void DumpyOS::buildIndex(DataSource* data_source) { // ---- searchIndex ---- +// Adapted from SaxUtil::getValueRange. +// Uses DaiSy's sax_breakpoints (same offset formula as Downloads' breakpoints array). +static void get_value_range(int sax_val, int bc, double* lb, double* ub) { + int cardinality = 1 << bc; + int offset = ((cardinality - 1) * (cardinality - 2)) / 2; + if (sax_val == 0) { + *lb = -std::numeric_limits::max(); + *ub = sax_breakpoints[offset]; + } else if (sax_val == cardinality - 1) { + *lb = sax_breakpoints[offset + sax_val - 1]; + *ub = std::numeric_limits::max(); + } else { + *lb = sax_breakpoints[offset + sax_val - 1]; + *ub = sax_breakpoints[offset + sax_val]; + } +} + +// Adapted from SaxUtil::LowerBound_Paa_iSax(paa, sax, bits_cardinality, chosen_segs, new_id). +// Computes the iSAX lower bound between paa and the potential child `child_id` of `node`. +static double lb_paa_to_child(const float* paa, const DumpyOSNode* node, + int child_id, int dim, int w) { + double coef = (double)dim / w; + double sum = 0.0; + int cid = child_id; + int cur = (int)node->chosen_segs.size() - 1; + + for (int i = w - 1; i >= 0; --i) { + int sax_val, bc; + if (cur >= 0 && node->chosen_segs[cur] == i) { + sax_val = (node->sax_word[i] << 1) | (cid & 1); + cid >>= 1; + bc = node->levels[i] + 1; + --cur; + } else { + sax_val = node->sax_word[i]; + bc = node->levels[i]; + } + if (bc == 0) continue; + double lo, hi; + get_value_range(sax_val, bc, &lo, &hi); + double p = paa[i]; + if (p < lo) sum += (lo - p) * (lo - p); + else if (p > hi) sum += (p - hi) * (p - hi); + } + return coef * sum; +} + static float l2sq_early(const float* a, const float* b, int d, float bound) { float s = 0.0f; for (int i = 0; i < d; ++i) { @@ -450,6 +507,7 @@ static float l2sq_early(const float* a, const float* b, int d, float bound) { return s; } +// Exact search adapted from FADASSearcher::exactSearchIdLevel. void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, idx_t* I, float* D) { if (!validateSearchParams(k, n_query)) return; @@ -462,51 +520,96 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, std::vector q_sax(w); std::vector q_paa(w); + struct PqItem { + double lb; + DumpyOSNode* parent; + int child_id; + bool operator>(const PqItem& o) const { return lb > o.lb; } + }; + for (idx_t qi = 0; qi < n_query; ++qi) { const float* q = query + (size_t)qi * dim; paa_from_ts(q, q_paa.data(), w, pts_per_seg); sax_from_paa(q_paa.data(), q_sax.data(), w, cardinality, max_bits); - // Route to leaf — equivalent to FADASNode::route() - DumpyOSNode* node = root_; - while (!node->chosen_segs.empty()) { - int sid = extend_sax_chosen(q_sax.data(), node->levels.data(), - node->chosen_segs, max_bits); - if (sid < 0 || sid >= (int)node->children.size() || - node->children[sid] == nullptr) break; - node = node->children[sid]; + // --- Approx phase: route to nearest leaf (FADASNode::route) --- + DumpyOSNode* approx_leaf = root_; + while (!approx_leaf->chosen_segs.empty()) { + int sid = extend_sax_chosen(q_sax.data(), approx_leaf->levels.data(), + approx_leaf->chosen_segs, max_bits); + if (sid < 0 || sid >= (int)approx_leaf->children.size() || + approx_leaf->children[sid] == nullptr) break; + approx_leaf = approx_leaf->children[sid]; } - // kNN via max-heap + // kNN max-heap: top = worst (largest) distance among current k-NN using Pair = std::pair; - std::vector heap; - heap.reserve((size_t)k + 1); + std::priority_queue> heap; float bsf = FLT_MAX; - for (idx_t si : node->entries) { - float dist = l2sq_early(q, database + (size_t)si * dim, (int)dim, bsf); - if ((idx_t)heap.size() < k || dist < bsf) { - heap.push_back({dist, si}); - std::push_heap(heap.begin(), heap.end()); - if ((idx_t)heap.size() > k) { - std::pop_heap(heap.begin(), heap.end()); - heap.pop_back(); + auto search_leaf = [&](DumpyOSNode* leaf) { + for (idx_t si : leaf->entries) { + float dist = l2sq_early(q, database + (size_t)si * dim, (int)dim, bsf); + if ((idx_t)heap.size() < k || dist < bsf) { + heap.push({dist, si}); + if ((idx_t)heap.size() > k) heap.pop(); + if ((idx_t)heap.size() == k) bsf = heap.top().first; } - if ((idx_t)heap.size() == k) bsf = heap.front().first; + } + }; + + if (approx_leaf->chosen_segs.empty()) + search_leaf(approx_leaf); + + // --- Exact phase (FADASSearcher::exactSearchIdLevel) --- + // Priority queue ordered ascending by iSAX lower bound. + std::priority_queue, std::greater> pq; + std::unordered_set visited; + visited.insert(approx_leaf); + + // Seed PQ with root's children (equivalent to seeding with vertexNum entries) + if (!root_->chosen_segs.empty()) { + for (int sid = 0; sid < (int)root_->children.size(); ++sid) { + DumpyOSNode* ch = root_->children[sid]; + if (ch == nullptr || visited.count(ch)) continue; + double lb = lb_paa_to_child(q_paa.data(), root_, sid, (int)dim, w); + pq.push({lb, root_, sid}); } } - std::sort_heap(heap.begin(), heap.end()); - for (idx_t j = 0; j < k; ++j) { - if (j < (idx_t)heap.size()) { - I[qi * k + j] = heap[j].second; - D[qi * k + j] = heap[j].first; + while (!pq.empty()) { + PqItem top = pq.top(); pq.pop(); + if (top.lb >= (double)bsf) break; + + DumpyOSNode* node = top.parent->children[top.child_id]; + if (node == nullptr || visited.count(node)) continue; + visited.insert(node); + + if (!node->chosen_segs.empty()) { + // Internal node: enqueue children with their LBs + for (int sid = 0; sid < (int)node->children.size(); ++sid) { + DumpyOSNode* ch = node->children[sid]; + if (ch == nullptr || visited.count(ch)) continue; + double lb = lb_paa_to_child(q_paa.data(), node, sid, (int)dim, w); + if (lb < (double)bsf) pq.push({lb, node, sid}); + } } else { - I[qi * k + j] = 0; - D[qi * k + j] = FLT_MAX; + search_leaf(node); } } + + // Extract results from max-heap in ascending order of distance + idx_t n_res = (idx_t)heap.size(); + for (idx_t j = n_res; j > 0; --j) { + I[qi * k + (j - 1)] = heap.top().second; + D[qi * k + (j - 1)] = heap.top().first; + heap.pop(); + } + for (idx_t j = n_res; j < k; ++j) { + I[qi * k + j] = 0; + D[qi * k + j] = FLT_MAX; + } } } diff --git a/lib/algos/DumpyOS.hpp b/lib/algos/DumpyOS.hpp index 8c003e8..f6f8c60 100644 --- a/lib/algos/DumpyOS.hpp +++ b/lib/algos/DumpyOS.hpp @@ -19,6 +19,7 @@ struct DumpyOSConfig { // Adapted from FADASNode in DumpyOS/include/DataStructures/FADASNode.h struct DumpyOSNode { std::vector levels; // bits_cardinality[] per segment + std::vector sax_word; // SAX word at current bit depth (needed for LB) std::vector chosen_segs; // chosenSegments; empty ↔ leaf std::vector children; // 2^|chosen_segs| entries (may be nullptr) std::vector entries; // series indices (leaf only) From dc870200b0218680f8f628e0bf16f99c1971d52c Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Fri, 26 Jun 2026 11:05:37 +0200 Subject: [PATCH 04/19] Add dumpy benchmark source and README --- README.md | 1 + benchmark/bm_DumpyOS_L2Square.cpp | 147 ++++++++++++++++++++++++++++++ 2 files changed, 148 insertions(+) create mode 100644 benchmark/bm_DumpyOS_L2Square.cpp diff --git a/README.md b/README.md index 935cab4..eaf84b6 100644 --- a/README.md +++ b/README.md @@ -46,6 +46,7 @@ The following table summarizes the key features of each algorithm: | **[Odyssey](https://helios2.mi.parisdescartes.fr/~themisp/odyssey/)** | Distributed and parallel in-memory similarity search | | **[SOFA](https://helios2.mi.parisdescartes.fr/~themisp/publications/icde25-sofa.pdf)** | In-memory similarity search using Symbolic Fourier Approximation (SFA) | | **[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 | diff --git a/benchmark/bm_DumpyOS_L2Square.cpp b/benchmark/bm_DumpyOS_L2Square.cpp new file mode 100644 index 0000000..2159743 --- /dev/null +++ b/benchmark/bm_DumpyOS_L2Square.cpp @@ -0,0 +1,147 @@ +#include +#include +#include +#include "bm_utils.hpp" +#include "../commons/dataloaders.hpp" +#include "../commons/VectorDataLoader.h" +#include "../commons/test_bm_utils.hpp" +#include "../lib/algos/DumpyOS.hpp" +#include "../lib/algos/DataSource.hpp" + +static bool endsWith(const std::string& s, const std::string& suffix) { + return s.size() >= suffix.size() && + s.compare(s.size() - suffix.size(), suffix.size(), suffix) == 0; +} + +struct DumpyOSSearchOnlyFixture : public benchmark::Fixture { + daisy::DumpyOS* search = nullptr; + float* database = nullptr; + float* query = nullptr; + daisy::idx_t* I = nullptr; + float* D = nullptr; + daisy::idx_t n_query = 0; + size_t k = 0; + std::string dataset_name; + size_t n_database = 0; + int thread_count = 0; + + void SetUp(const benchmark::State& state) override { + int config_idx = static_cast(state.range(0)); + const SSTestConfig& config = test_configs_deep_seismic_astro270m[config_idx]; + + const bool use_fvecs = endsWith(config.dataset_path, ".fvecs") || endsWith(config.query_path, ".fvecs"); + size_t dim_u = 0, n_database_u = 0, n_q_u = 0; + database = nullptr; + + if (use_fvecs) { + database = fvecs_read(config.dataset_path.c_str(), &dim_u, &n_database_u, 0); + if (!database) { + std::cerr << "Failed to load dataset (fvecs)" << std::endl; + return; + } + const size_t query_limit = (config.query_limit > 0) ? static_cast(config.query_limit) : 0; + query = fvecs_read(config.query_path.c_str(), &dim_u, &n_q_u, query_limit); + if (!query) { + std::cerr << "Failed to load queries (fvecs)" << std::endl; + delete[] database; + return; + } + } else { + std::string dataset_filename = pathToFilename(config.dataset_path); + std::string query_filename = pathToFilename(config.query_path); + + daisy::idx_t dim, n_database, _, __; + if (!parseFilenameForConfig(dataset_filename, "bruteForce", dim, n_database, _, __)) { + std::cerr << "Failed to parse dataset config from filename: " << dataset_filename << std::endl; + return; + } + + daisy::idx_t dim_q, n_q, ___, ____; + if (!parseFilenameForConfig(query_filename, "bruteForce", dim_q, n_q, ___, ____)) { + std::cerr << "Failed to parse query config from filename: " << query_filename << std::endl; + return; + } + + if (dim != static_cast(dim_q)) { + std::cerr << "Dimension mismatch between dataset and queries" << std::endl; + return; + } + + dim_u = static_cast(dim); + n_database_u = static_cast(n_database); + if (config.query_limit > 0 && static_cast(config.query_limit) < n_q) + n_q = static_cast(config.query_limit); + n_q_u = static_cast(n_q); + + database = loadBinData(config.dataset_path.c_str(), n_database, dim, false); + if (!database) { + std::cerr << "Failed to load dataset" << std::endl; + return; + } + + query = loadBinData(config.query_path.c_str(), n_q, dim_q, false); + if (!query) { + std::cerr << "Failed to load queries" << std::endl; + delete[] database; + return; + } + } + + search = new daisy::DumpyOS(daisy::DistanceType::L2_SQUARED); + search->setNumThreads(config.thread_count); + + fprintf(stderr, "[DUMPYOS] Before buildIndex (n_database=%zu dim=%zu).\n", n_database_u, dim_u); + fflush(stderr); + + daisy::InMemoryDataSource data_source(database, static_cast(n_database_u), static_cast(dim_u)); + search->buildIndex(&data_source); + + fprintf(stderr, "[DUMPYOS] Indexing finished (n_database=%zu dim=%zu).\n", n_database_u, dim_u); + fflush(stderr); + + k = static_cast(config.k_value); + n_query = static_cast(n_q_u); + I = new daisy::idx_t[n_query * k]; + D = new float[n_query * k]; + + dataset_name = config.name; + n_database = n_database_u; + thread_count = config.thread_count; + + fprintf(stderr, "[DUMPYOS] n_database=%zu n_query=%zu dim=%zu k=%zu threads=%d\n", + n_database_u, (size_t)n_query, dim_u, k, config.thread_count); + fflush(stderr); + } + + void TearDown(const benchmark::State&) override { + delete search; + delete[] database; + delete[] query; + delete[] I; + delete[] D; + search = nullptr; + database = nullptr; + query = nullptr; + I = nullptr; + D = nullptr; + } +}; + +BENCHMARK_DEFINE_F(DumpyOSSearchOnlyFixture, BM_DumpyOS_SearchOnly)(benchmark::State& state) { + for (auto _ : state) { + fprintf(stderr, "[DUMPYOS] --- Query phase ---\n"); + fprintf(stderr, "[DUMPYOS] dataset=%s n_database=%zu\n", dataset_name.c_str(), n_database); + fprintf(stderr, "[DUMPYOS] search_threads=%d n_query=%zu k=%zu\n", thread_count, (size_t)n_query, k); + fflush(stderr); + search->searchIndex(query, n_query, static_cast(k), I, D); + fprintf(stderr, "[DUMPYOS] Querying finished (n_query=%zu k=%zu).\n", (size_t)n_query, k); + fflush(stderr); + } +} + +BENCHMARK_REGISTER_F(DumpyOSSearchOnlyFixture, BM_DumpyOS_SearchOnly) + ->Args({0})->Args({1})->Args({2})->Args({3})->Args({4})->Args({5})->Args({6})->Args({7}) + ->Iterations(1) + ->Unit(benchmark::kMillisecond); + +BENCHMARK_MAIN(); From 8a6244d7006a5aea5df3836a096d028bf1ae6b42 Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Fri, 26 Jun 2026 11:17:32 +0200 Subject: [PATCH 05/19] Final dumpy (with passed tests) --- lib/algos/DumpyOS.cpp | 88 +++++++++---------------------------------- 1 file changed, 17 insertions(+), 71 deletions(-) diff --git a/lib/algos/DumpyOS.cpp b/lib/algos/DumpyOS.cpp index 5e740e9..2dd0984 100644 --- a/lib/algos/DumpyOS.cpp +++ b/lib/algos/DumpyOS.cpp @@ -1,8 +1,5 @@ -// Adapted from FADASNode.cpp in DumpyOS/src/IndexConstruction/FADASNode.cpp -// (Wang et al., "DumpyOS: A Dumpy Index for Scalable Data Series Similarity -// Search", VLDB Journal 2024) -// Reference: SaxUtil.cpp (getMidLineFromSaxSymbolbc8, extendSax), -// MathUtil.cpp (deviation, generateMaskSettingKbits) +// based on FADASNode.cpp, SaxUtil.cpp and MathUtil.cpp from the DumpyOS repo +// (Wang et al., "DumpyOS: A Dumpy Index for Scalable Data Series Similarity Search", VLDB Journal 2024) #include "DumpyOS.hpp" #include "../isax/SAX.hpp" @@ -19,10 +16,7 @@ namespace daisy { -// ---- bp8 : normal-distribution breakpoints for 8-bit SAX ---- -// Directly from SaxUtil.cpp in the DumpyOS reference implementation. -// bp8[i] is the upper boundary of the i-th SAX symbol region. -// bp8[255] = +inf sentinel; midpoints are extrapolated for symbols 0 and 255. +// normal distribution breakpoints for 8-bit SAX (from SaxUtil.cpp) static const double bp8[256] = { -2.660067468617458,-2.4175590162365035,-2.2662268092096522, -2.1538746940614573,-2.063527898316245,-1.9874278859298962, @@ -111,21 +105,17 @@ static const double bp8[256] = { 1.921350774293703, 1.9874278859298962, 2.063527898316245, 2.1538746940614573, 2.2662268092096522, 2.4175590162365035, 2.660067468617458, - // sentinel for symbol 255 upper bound (numeric_limits::max()) 3.4028234663852886e+38 }; -// ---- internal helpers (adapted from SaxUtil and MathUtil) ---- - -// Equivalent to SaxUtil::getMidLineFromSaxSymbolbc8. +// midpoint of the SAX symbol region (getMidLineFromSaxSymbolbc8) static double get_midpoint(unsigned sym) { if (sym == 0) return bp8[0] - (bp8[1] - bp8[0]); if (sym == 255) return bp8[254] + (bp8[254] - bp8[253]); return (bp8[sym - 1] + bp8[sym]) * 0.5; } -// Equivalent to SaxUtil::extendSax(sax, bits_cardinality) — all segments. -// Produces a vertexNum (1<& chosen, int max_bits) { int res = 0; @@ -149,8 +138,7 @@ static int extend_sax_chosen(const sax_type* sax, const int* levels, return res; } -// Equivalent to MathUtil::generateMaskSettingKbits(plan, k, len). -// plan must be sorted ascending. Sets bit (w-1-plan[i]) for each i in plan. +// bitmask with 1s at positions corresponding to the chosen segments static int generate_mask(const int* plan, int lambda, int w) { int res = 0, cur = 0; for (int i = 0; i < w; ++i) { @@ -160,7 +148,6 @@ static int generate_mask(const int* plan, int lambda, int w) { return res; } -// Population standard deviation — equivalent to MathUtil::deviation(double*, int). static double pop_stdev(const double* vals, int n) { double mean = 0.0; for (int i = 0; i < n; ++i) mean += vals[i]; @@ -170,9 +157,6 @@ static double pop_stdev(const double* vals, int n) { return std::sqrt(var / n); } -// Equivalent to FADASNode::compute_score. -// node_n = node->n (total series in node) -// alpha = config_.alpha static double compute_score_(const std::vector& node_sizes, const int* plan, int lambda, const std::vector& data_seg_stdev, @@ -199,10 +183,7 @@ static double compute_score_(const std::vector& node_sizes, return sum_seg + alpha * balance; } -// Equivalent to FADASNode::visitPlanFromBaseTable. -// cur_lambda = parent lambda minus 1 (we're evaluating plans of size cur_lambda). -// plan = parent plan of size cur_lambda + 1 (sorted ascending). -// base_tbl = parent child fill counts of size 1 << (cur_lambda + 1). +// recursively evaluates sub-plans by dropping one segment at a time from the parent plan static void visit_plan_from_base_table( std::unordered_set& visited, int cur_lambda, const int* plan, const std::vector& base_tbl, @@ -211,22 +192,19 @@ static void visit_plan_from_base_table( const std::vector& data_seg_stdev, int node_n, int leaf_size, double alpha, int w) { - // base_mask has (cur_lambda+1) bits set (covers all 1<<(cur_lambda+1) entries) int base_mask = 1; for (int i = 0; i < cur_lambda; ++i) base_mask = (base_mask << 1) | 1; for (int i = 0; i <= cur_lambda; ++i) { - int reset_pos = plan[i]; - int cur_whole_mask = mask_code - (1 << (w - 1 - reset_pos)); + int reset_pos = plan[i]; + int cur_whole_mask = mask_code - (1 << (w - 1 - reset_pos)); if (visited.count(cur_whole_mask)) continue; visited.insert(cur_whole_mask); - // Sub-plan: drop element i from parent plan int* new_plan = new int[cur_lambda]; for (int j = 0, k = 0; j <= cur_lambda; ++j) if (j != i) new_plan[k++] = plan[j]; - // Aggregate base_tbl into cur_lambda child fill counts by masking out bit i int cur_base_mask = base_mask - (1 << (cur_lambda - i)); std::map nsmap; for (int j = 0; j < (int)base_tbl.size(); ++j) @@ -250,7 +228,6 @@ static void visit_plan_from_base_table( } } -// Equivalent to FADASNode::determineFanout. static void determine_fanout(int n, int leaf_size, double f_low, double f_high, int w, int* lambda_min, int* lambda_max) { if (n < 2 * leaf_size) { *lambda_min = 1; *lambda_max = 1; return; } @@ -271,7 +248,6 @@ static void determine_fanout(int n, int leaf_size, double f_low, double f_high, if (*lambda_min == -1) { *lambda_min = w; *lambda_max = w; } } -// Equivalent to FADASNode::determineSegments. void DumpyOS::determineSegments_(DumpyOSNode* node) { int w = config_.paa_segments; int max_bits = config_.sax_bit_cardinality; @@ -282,7 +258,6 @@ void DumpyOS::determineSegments_(DumpyOSNode* node) { determine_fanout(node->n, ls, config_.fill_lower, config_.fill_upper, w, &lambda_min, &lambda_max); - // Edge case: need all segments (vertexNum < _min) if (lambda_min == w && lambda_max == w) { for (int i = 0; i < w; ++i) node->chosen_segs.push_back(i); return; @@ -290,9 +265,7 @@ void DumpyOS::determineSegments_(DumpyOSNode* node) { int vertex_num = 1 << w; - // Build unit_size table: aggregate series counts into 2^w first-level buckets std::vector unit_size(vertex_num, 0); - // Per-segment SAX symbol frequency (for data_seg_stdev) std::vector> seg_sym_cnt(w); for (idx_t si : node->entries) { @@ -303,8 +276,6 @@ void DumpyOS::determineSegments_(DumpyOSNode* node) { unit_size[head]++; } - // Compute data_seg_mean and data_seg_stdev (variance of midpoint values) - // Equivalent to the data_seg_mean/stdev loop in FADASNode::determineSegments. std::vector data_seg_mean(w, 0.0), data_seg_stdev(w, 0.0); for (int i = 0; i < w; ++i) { for (auto& kv : seg_sym_cnt[i]) @@ -317,19 +288,16 @@ void DumpyOS::determineSegments_(DumpyOSNode* node) { data_seg_stdev[i] /= node->n; } - double max_score = 0.0; - std::vector best_plan; + double max_score = 0.0; + std::vector best_plan; std::unordered_set visited; - // Enumerate all C(w, lambda_max) plans; for each, also search sub-plans - // via visitPlanFromBaseTable (lambda_min..lambda_max-1). std::vector idx(lambda_max); std::iota(idx.begin(), idx.end(), 0); while (true) { const int* plan = idx.data(); - // Compute child fill counts using unit_size table + bit masking int mask_code = generate_mask(plan, lambda_max, w); std::map node_size_map; for (int j = 0; j < vertex_num; ++j) @@ -350,7 +318,7 @@ void DumpyOS::determineSegments_(DumpyOSNode* node) { &max_score, best_plan, lambda_min, mask_code, data_seg_stdev, node->n, ls, alpha, w); - // Advance to next C(w, lambda_max) combination (Gosper-style) + // next combination (Gosper-style) int i = lambda_max - 1; while (i >= 0 && idx[i] == w - lambda_max + i) --i; if (i < 0) break; @@ -361,26 +329,22 @@ void DumpyOS::determineSegments_(DumpyOSNode* node) { node->chosen_segs = best_plan; } -// ---- splitNode_ ---- - void DumpyOS::splitNode_(DumpyOSNode* node) { int w = config_.paa_segments; int max_bits = config_.sax_bit_cardinality; determineSegments_(node); - if (node->chosen_segs.empty()) return; // exhausted bit budget + if (node->chosen_segs.empty()) return; int lam = (int)node->chosen_segs.size(); int num_ch = 1 << lam; node->children.assign(num_ch, nullptr); - // Create all children with inherited levels/sax_word (+ 1 bit for each chosen segment) for (int sid = 0; sid < num_ch; ++sid) { DumpyOSNode* child = new DumpyOSNode(); child->levels = node->levels; child->sax_word = node->sax_word; int cid = sid; - // Extract bits LSB-first: last chosen_seg gets bit 0 of sid for (int j = lam - 1; j >= 0; --j) { int cs = node->chosen_segs[j]; child->sax_word[cs] = (child->sax_word[cs] << 1) | (cid & 1); @@ -390,7 +354,6 @@ void DumpyOS::splitNode_(DumpyOSNode* node) { node->children[sid] = child; } - // Redistribute series — equivalent to FADASNode's child routing for (idx_t si : node->entries) { const sax_type* sax = sax_table_ + (size_t)si * w; int sid = extend_sax_chosen(sax, node->levels.data(), node->chosen_segs, max_bits); @@ -405,8 +368,6 @@ void DumpyOS::splitNode_(DumpyOSNode* node) { splitNode_(child); } -// ---- buildIndex ---- - void DumpyOS::buildIndex(DataSource* data_source) { dim = data_source->getDim(); n_database = data_source->getTotalRecords(); @@ -427,7 +388,6 @@ void DumpyOS::buildIndex(DataSource* data_source) { int cardinality = 1 << max_bits; int pts_per_seg = (int)dim / w; - // Compute and cache the SAX word for every database series sax_table_ = new sax_type[(size_t)n_database * w]; std::vector paa(w); @@ -436,7 +396,6 @@ void DumpyOS::buildIndex(DataSource* data_source) { sax_from_paa(paa.data(), sax_table_ + (size_t)i * w, w, cardinality, max_bits); } - // Root: all series, all segment levels and sax_word at 0 root_ = new DumpyOSNode(); root_->levels.assign(w, 0); root_->sax_word.assign(w, 0); @@ -448,10 +407,7 @@ void DumpyOS::buildIndex(DataSource* data_source) { splitNode_(root_); } -// ---- searchIndex ---- - -// Adapted from SaxUtil::getValueRange. -// Uses DaiSy's sax_breakpoints (same offset formula as Downloads' breakpoints array). +// iSAX lower bound for a potential child node (LowerBound_Paa_iSax) static void get_value_range(int sax_val, int bc, double* lb, double* ub) { int cardinality = 1 << bc; int offset = ((cardinality - 1) * (cardinality - 2)) / 2; @@ -467,8 +423,6 @@ static void get_value_range(int sax_val, int bc, double* lb, double* ub) { } } -// Adapted from SaxUtil::LowerBound_Paa_iSax(paa, sax, bits_cardinality, chosen_segs, new_id). -// Computes the iSAX lower bound between paa and the potential child `child_id` of `node`. static double lb_paa_to_child(const float* paa, const DumpyOSNode* node, int child_id, int dim, int w) { double coef = (double)dim / w; @@ -507,7 +461,6 @@ static float l2sq_early(const float* a, const float* b, int d, float bound) { return s; } -// Exact search adapted from FADASSearcher::exactSearchIdLevel. void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, idx_t* I, float* D) { if (!validateSearchParams(k, n_query)) return; @@ -533,7 +486,7 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, paa_from_ts(q, q_paa.data(), w, pts_per_seg); sax_from_paa(q_paa.data(), q_sax.data(), w, cardinality, max_bits); - // --- Approx phase: route to nearest leaf (FADASNode::route) --- + // route to nearest leaf DumpyOSNode* approx_leaf = root_; while (!approx_leaf->chosen_segs.empty()) { int sid = extend_sax_chosen(q_sax.data(), approx_leaf->levels.data(), @@ -543,7 +496,6 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, approx_leaf = approx_leaf->children[sid]; } - // kNN max-heap: top = worst (largest) distance among current k-NN using Pair = std::pair; std::priority_queue> heap; float bsf = FLT_MAX; @@ -562,13 +514,11 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, if (approx_leaf->chosen_segs.empty()) search_leaf(approx_leaf); - // --- Exact phase (FADASSearcher::exactSearchIdLevel) --- - // Priority queue ordered ascending by iSAX lower bound. + // exact phase: priority queue pruned by iSAX lower bounds std::priority_queue, std::greater> pq; std::unordered_set visited; visited.insert(approx_leaf); - // Seed PQ with root's children (equivalent to seeding with vertexNum entries) if (!root_->chosen_segs.empty()) { for (int sid = 0; sid < (int)root_->children.size(); ++sid) { DumpyOSNode* ch = root_->children[sid]; @@ -587,7 +537,6 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, visited.insert(node); if (!node->chosen_segs.empty()) { - // Internal node: enqueue children with their LBs for (int sid = 0; sid < (int)node->children.size(); ++sid) { DumpyOSNode* ch = node->children[sid]; if (ch == nullptr || visited.count(ch)) continue; @@ -599,7 +548,6 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, } } - // Extract results from max-heap in ascending order of distance idx_t n_res = (idx_t)heap.size(); for (idx_t j = n_res; j > 0; --j) { I[qi * k + (j - 1)] = heap.top().second; @@ -613,8 +561,6 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, } } -// ---- destructor ---- - void DumpyOS::destroyTree_(DumpyOSNode* node) { if (!node) return; for (DumpyOSNode* child : node->children) destroyTree_(child); From cc59c1f93e14cf5ab56cde0f5c4bdd634ec17c9b Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Fri, 26 Jun 2026 14:55:18 +0200 Subject: [PATCH 06/19] DTW support for DumpyOS --- lib/algos/DumpyOS.cpp | 73 +++++++++++++++++++++++++++++++++----- tests/CMakeLists.txt | 27 ++++++++++++++ tests/test_DumpyOS_DTW.cpp | 41 +++++++++++++++++++++ tests/test_utils.hpp | 10 ++++++ 4 files changed, 143 insertions(+), 8 deletions(-) create mode 100644 tests/test_DumpyOS_DTW.cpp diff --git a/lib/algos/DumpyOS.cpp b/lib/algos/DumpyOS.cpp index 2dd0984..30260df 100644 --- a/lib/algos/DumpyOS.cpp +++ b/lib/algos/DumpyOS.cpp @@ -1,8 +1,8 @@ -// based on FADASNode.cpp, SaxUtil.cpp and MathUtil.cpp from the DumpyOS repo -// (Wang et al., "DumpyOS: A Dumpy Index for Scalable Data Series Similarity Search", VLDB Journal 2024) #include "DumpyOS.hpp" #include "../isax/SAX.hpp" +#include "../isax/iSAXIndex.hpp" +#include "../distance_computers/DistanceComputer.hpp" #include #include @@ -16,7 +16,7 @@ namespace daisy { -// normal distribution breakpoints for 8-bit SAX (from SaxUtil.cpp) +// normal distribution breakpoints for 8-bit SAX (from SaxUtil.cpp, og DUMPYOS code) static const double bp8[256] = { -2.660067468617458,-2.4175590162365035,-2.2662268092096522, -2.1538746940614573,-2.063527898316245,-1.9874278859298962, @@ -183,7 +183,7 @@ static double compute_score_(const std::vector& node_sizes, return sum_seg + alpha * balance; } -// recursively evaluates sub-plans by dropping one segment at a time from the parent plan +// recursively evaluate sub-plans by dropping one segment at a time from the parent plan static void visit_plan_from_base_table( std::unordered_set& visited, int cur_lambda, const int* plan, const std::vector& base_tbl, @@ -451,6 +451,35 @@ static double lb_paa_to_child(const float* paa, const DumpyOSNode* node, return coef * sum; } +// DTW version of lb_paa_to_child: uses [paaL, paaU] interval instead of single paa value (again, from dumpyos implementation) +static double lb_paa_to_child_dtw(const float* paaU, const float* paaL, + const DumpyOSNode* node, + int child_id, int dim, int w) { + double coef = (double)dim / w; + double sum = 0.0; + int cid = child_id; + int cur = (int)node->chosen_segs.size() - 1; + + for (int i = w - 1; i >= 0; --i) { + int sax_val, bc; + if (cur >= 0 && node->chosen_segs[cur] == i) { + sax_val = (node->sax_word[i] << 1) | (cid & 1); + cid >>= 1; + bc = node->levels[i] + 1; + --cur; + } else { + sax_val = node->sax_word[i]; + bc = node->levels[i]; + } + if (bc == 0) continue; + double lo, hi; + get_value_range(sax_val, bc, &lo, &hi); + if (paaU[i] < lo) sum += (lo - paaU[i]) * (lo - paaU[i]); + else if (paaL[i] > hi) sum += (paaL[i] - hi) * (paaL[i] - hi); + } + return coef * sum; +} + static float l2sq_early(const float* a, const float* b, int d, float bound) { float s = 0.0f; for (int i = 0; i < d; ++i) { @@ -465,13 +494,19 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, idx_t* I, float* D) { if (!validateSearchParams(k, n_query)) return; + bool use_dtw = (this->distance_type == DistanceType::DTW); int w = config_.paa_segments; int max_bits = config_.sax_bit_cardinality; int cardinality = 1 << max_bits; int pts_per_seg = (int)dim / w; + int warp_win = static_cast(dim * 0.1); std::vector q_sax(w); std::vector q_paa(w); + std::vector q_paa_upper(use_dtw ? w : 0); + std::vector q_paa_lower(use_dtw ? w : 0); + std::vector upper_env(use_dtw ? (int)dim : 0); + std::vector lower_env(use_dtw ? (int)dim : 0); struct PqItem { double lb; @@ -486,6 +521,13 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, paa_from_ts(q, q_paa.data(), w, pts_per_seg); sax_from_paa(q_paa.data(), q_sax.data(), w, cardinality, max_bits); + if (use_dtw) { + lower_upper_lemire(const_cast(q), (int)dim, warp_win, + lower_env.data(), upper_env.data()); + paa_from_ts(upper_env.data(), q_paa_upper.data(), w, pts_per_seg); + paa_from_ts(lower_env.data(), q_paa_lower.data(), w, pts_per_seg); + } + // route to nearest leaf DumpyOSNode* approx_leaf = root_; while (!approx_leaf->chosen_segs.empty()) { @@ -502,7 +544,15 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, auto search_leaf = [&](DumpyOSNode* leaf) { for (idx_t si : leaf->entries) { - float dist = l2sq_early(q, database + (size_t)si * dim, (int)dim, bsf); + float dist; + if (use_dtw) { + dist = distance_computer->compute_dist( + const_cast(q), + database + (size_t)si * dim, + (int)dim, bsf); + } else { + dist = l2sq_early(q, database + (size_t)si * dim, (int)dim, bsf); + } if ((idx_t)heap.size() < k || dist < bsf) { heap.push({dist, si}); if ((idx_t)heap.size() > k) heap.pop(); @@ -511,6 +561,13 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, } }; + auto compute_lb = [&](const DumpyOSNode* parent, int sid) -> double { + if (use_dtw) + return lb_paa_to_child_dtw(q_paa_upper.data(), q_paa_lower.data(), + parent, sid, (int)dim, w); + return lb_paa_to_child(q_paa.data(), parent, sid, (int)dim, w); + }; + if (approx_leaf->chosen_segs.empty()) search_leaf(approx_leaf); @@ -523,7 +580,7 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, for (int sid = 0; sid < (int)root_->children.size(); ++sid) { DumpyOSNode* ch = root_->children[sid]; if (ch == nullptr || visited.count(ch)) continue; - double lb = lb_paa_to_child(q_paa.data(), root_, sid, (int)dim, w); + double lb = compute_lb(root_, sid); pq.push({lb, root_, sid}); } } @@ -540,7 +597,7 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, for (int sid = 0; sid < (int)node->children.size(); ++sid) { DumpyOSNode* ch = node->children[sid]; if (ch == nullptr || visited.count(ch)) continue; - double lb = lb_paa_to_child(q_paa.data(), node, sid, (int)dim, w); + double lb = compute_lb(node, sid); if (lb < (double)bsf) pq.push({lb, node, sid}); } } else { @@ -579,4 +636,4 @@ DumpyOS::~DumpyOS() { if (owns_database_) delete[] database; } -} // namespace daisy +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 10eb86a..8ad95bd 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -651,4 +651,31 @@ target_include_directories(test_DumpyOS_L2Square gtest_discover_tests( test_DumpyOS_L2Square WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} +) + +# ////// DUMPYOS DTW ////// +add_executable( + test_DumpyOS_DTW + test_DumpyOS_DTW.cpp + test_utils.cpp +) + +target_link_libraries( + test_DumpyOS_DTW + PRIVATE + GTest::gtest_main + dino_lib + commons_lib + stdc++fs +) + +target_include_directories(test_DumpyOS_DTW + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons +) + +gtest_discover_tests( + test_DumpyOS_DTW + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} ) \ No newline at end of file diff --git a/tests/test_DumpyOS_DTW.cpp b/tests/test_DumpyOS_DTW.cpp new file mode 100644 index 0000000..50dbd67 --- /dev/null +++ b/tests/test_DumpyOS_DTW.cpp @@ -0,0 +1,41 @@ +#include "test_utils.hpp" +#include "../commons/test_bm_utils.hpp" +#include "../commons/paramSetup.hpp" + +std::string prefix = "bruteForce"; + +TEST_P(DumpyOSDTWParameterizedTest, AllConfigurations) +{ + const SSTestConfig &config = GetParam(); + daisy::DistanceType dist_DTW = daisy::DistanceType::DTW; + daisy::DumpyOS search(dist_DTW); + + std::string gt_I_path = config.gt_I_prefix + std::to_string(config.k_value) + ".txt"; + std::string gt_D_path = config.gt_D_prefix + std::to_string(config.k_value) + ".txt"; + + runSST( + &search, + prefix, + gt_I_path, + gt_D_path, + config.dataset_path, + config.query_path, + config.thread_count); +} + +INSTANTIATE_TEST_SUITE_P( + DumpyOSDTWTests, + DumpyOSDTWParameterizedTest, + ::testing::ValuesIn(test_configs_dtw), + [](const ::testing::TestParamInfo &info) + { + return info.param.name + "_k" + std::to_string(info.param.k_value) + + "_thread" + std::to_string(info.param.thread_count) + + "_idx" + std::to_string(info.index); + }); + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/tests/test_utils.hpp b/tests/test_utils.hpp index 3b52b2b..8ac8a56 100644 --- a/tests/test_utils.hpp +++ b/tests/test_utils.hpp @@ -153,4 +153,14 @@ class DumpyOSParameterizedTest : public SimilaritySearchTest, static void TearDownTestSuite() {} }; +class DumpyOSDTWParameterizedTest : public SimilaritySearchTest, + public ::testing::WithParamInterface +{ +protected: + using SimilaritySearchTest::runSST; + + static void SetUpTestSuite() {} + static void TearDownTestSuite() {} +}; + #endif From 63d3ecfc3522735e1f69f6513752a54cb1aca13a Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Fri, 26 Jun 2026 15:30:20 +0200 Subject: [PATCH 07/19] fix DTW bug --- lib/algos/DumpyOS.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/lib/algos/DumpyOS.cpp b/lib/algos/DumpyOS.cpp index 30260df..3b894d7 100644 --- a/lib/algos/DumpyOS.cpp +++ b/lib/algos/DumpyOS.cpp @@ -2,6 +2,7 @@ #include "DumpyOS.hpp" #include "../isax/SAX.hpp" #include "../isax/iSAXIndex.hpp" +#include "../isax/iSAXSearch.hpp" #include "../distance_computers/DistanceComputer.hpp" #include @@ -507,6 +508,8 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, std::vector q_paa_lower(use_dtw ? w : 0); std::vector upper_env(use_dtw ? (int)dim : 0); std::vector lower_env(use_dtw ? (int)dim : 0); + // all-zeros cb: dtw() only reads cb[i+r+1] for pruning; with zeros it's equivalent to no pruning + std::vector dtw_cb(use_dtw ? (int)dim : 0, 0.0f); struct PqItem { double lb; @@ -546,10 +549,10 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, for (idx_t si : leaf->entries) { float dist; if (use_dtw) { - dist = distance_computer->compute_dist( - const_cast(q), - database + (size_t)si * dim, - (int)dim, bsf); + // use naive dtw (correct), dtwsimdPruned has a buggy DP initialization + dist = dtw(const_cast(q), + database + (size_t)si * dim, + dtw_cb.data(), (int)dim, warp_win, bsf); } else { dist = l2sq_early(q, database + (size_t)si * dim, (int)dim, bsf); } From 6a1bdec33bebad05a8aa808f3b8ac8090f8af60b Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Fri, 26 Jun 2026 15:46:38 +0200 Subject: [PATCH 08/19] revert dtw leaf distance to compute_dist (same as rest of daisy) --- lib/algos/DumpyOS.cpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/lib/algos/DumpyOS.cpp b/lib/algos/DumpyOS.cpp index 3b894d7..30260df 100644 --- a/lib/algos/DumpyOS.cpp +++ b/lib/algos/DumpyOS.cpp @@ -2,7 +2,6 @@ #include "DumpyOS.hpp" #include "../isax/SAX.hpp" #include "../isax/iSAXIndex.hpp" -#include "../isax/iSAXSearch.hpp" #include "../distance_computers/DistanceComputer.hpp" #include @@ -508,8 +507,6 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, std::vector q_paa_lower(use_dtw ? w : 0); std::vector upper_env(use_dtw ? (int)dim : 0); std::vector lower_env(use_dtw ? (int)dim : 0); - // all-zeros cb: dtw() only reads cb[i+r+1] for pruning; with zeros it's equivalent to no pruning - std::vector dtw_cb(use_dtw ? (int)dim : 0, 0.0f); struct PqItem { double lb; @@ -549,10 +546,10 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, for (idx_t si : leaf->entries) { float dist; if (use_dtw) { - // use naive dtw (correct), dtwsimdPruned has a buggy DP initialization - dist = dtw(const_cast(q), - database + (size_t)si * dim, - dtw_cb.data(), (int)dim, warp_win, bsf); + dist = distance_computer->compute_dist( + const_cast(q), + database + (size_t)si * dim, + (int)dim, bsf); } else { dist = l2sq_early(q, database + (size_t)si * dim, (int)dim, bsf); } From 59a4b0ac2f241e39e6c9371b069cd75a02a8efd2 Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Tue, 30 Jun 2026 10:42:14 +0200 Subject: [PATCH 09/19] add demo for dumpy dtw --- demos/CMakeLists.txt | 30 ++++++++++++++++++++++++++ demos/demo_DumpyOS_DTW.cpp | 44 ++++++++++++++++++++++++++++++++++++++ demos/demo_DumpyOS_DTW.py | 37 ++++++++++++++++++++++++++++++++ lib/algos/DumpyOS.cpp | 26 +++++++++++----------- lib/algos/DumpyOS.hpp | 8 ++++--- pybinds/setup.cpp | 1 + 6 files changed, 130 insertions(+), 16 deletions(-) create mode 100644 demos/demo_DumpyOS_DTW.cpp create mode 100644 demos/demo_DumpyOS_DTW.py diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt index 6e330d7..f85c069 100644 --- a/demos/CMakeLists.txt +++ b/demos/CMakeLists.txt @@ -466,6 +466,36 @@ if(BUILD_DEMO) message(STATUS "Include directories added for demo_DumpyOS_L2Square.") endif() + # ////// DUMPYOS DTW ////// + if(DEBUG_MSG) + message(STATUS "---") + message(STATUS "## Demo: DumpyOS DTW") + message(STATUS "Attempting to add executable: demo_DumpyOS_DTW") + endif() + add_executable(demo_DumpyOS_DTW demo_DumpyOS_DTW.cpp) + if(DEBUG_MSG) + message(STATUS "Executable demo_DumpyOS_DTW added.") + endif() + + if(DEBUG_MSG) + message(STATUS "Linking libraries for demo_DumpyOS_DTW...") + endif() + target_link_libraries(demo_DumpyOS_DTW PRIVATE dino_lib commons_lib) + if(DEBUG_MSG) + message(STATUS "Libraries linked for demo_DumpyOS_DTW.") + endif() + + if(DEBUG_MSG) + message(STATUS "Adding include directories for demo_DumpyOS_DTW...") + endif() + target_include_directories(demo_DumpyOS_DTW PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + if(DEBUG_MSG) + message(STATUS "Include directories added for demo_DumpyOS_DTW.") + endif() + else() if(DEBUG_MSG) message(STATUS "BUILD_DEMO is FALSE. Demo executables will NOT be built.") diff --git a/demos/demo_DumpyOS_DTW.cpp b/demos/demo_DumpyOS_DTW.cpp new file mode 100644 index 0000000..61565f0 --- /dev/null +++ b/demos/demo_DumpyOS_DTW.cpp @@ -0,0 +1,44 @@ +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" +#include +#include + +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::DumpyOS dumpyos_search(daisy::DistanceType::DTW); + dumpyos_search.setNumThreads(4); + + int warp_window = std::max(1, static_cast(dim * 0.1)); + dumpyos_search.setWarpingWindow(warp_window); + + dumpyos_search.buildIndex(database, n_database, dim); + + daisy::idx_t *I = new daisy::idx_t[n_query * k]; + float *D = new float[n_query * k]; + dumpyos_search.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_DumpyOS_DTW.py b/demos/demo_DumpyOS_DTW.py new file mode 100644 index 0000000..0ac8bc9 --- /dev/null +++ b/demos/demo_DumpyOS_DTW.py @@ -0,0 +1,37 @@ +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, DumpyOS + +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 = DumpyOS(DistanceType.DTW) + index.setWarpingWindow(max(1, int(dim * 0.1))) + + index.setNumThreads(4) + index.buildIndex(db) + + 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/lib/algos/DumpyOS.cpp b/lib/algos/DumpyOS.cpp index 30260df..f87bdd9 100644 --- a/lib/algos/DumpyOS.cpp +++ b/lib/algos/DumpyOS.cpp @@ -494,24 +494,24 @@ void DumpyOS::searchIndex(const float* query, idx_t n_query, idx_t k, idx_t* I, float* D) { if (!validateSearchParams(k, n_query)) return; - bool use_dtw = (this->distance_type == DistanceType::DTW); - int w = config_.paa_segments; - int max_bits = config_.sax_bit_cardinality; + bool use_dtw = (this->distance_type == DistanceType::DTW); + int w = config_.paa_segments; + int max_bits = config_.sax_bit_cardinality; int cardinality = 1 << max_bits; int pts_per_seg = (int)dim / w; - int warp_win = static_cast(dim * 0.1); + int warp_win = warping_window; std::vector q_sax(w); - std::vector q_paa(w); - std::vector q_paa_upper(use_dtw ? w : 0); - std::vector q_paa_lower(use_dtw ? w : 0); - std::vector upper_env(use_dtw ? (int)dim : 0); - std::vector lower_env(use_dtw ? (int)dim : 0); - - struct PqItem { - double lb; + std::vector q_paa(w); + std::vector q_paa_upper(use_dtw ? w : 0); + std::vector q_paa_lower(use_dtw ? w : 0); + std::vector upper_env(use_dtw ? (int)dim : 0); + std::vector lower_env(use_dtw ? (int)dim : 0); + + struct PqItem{ + double lb; DumpyOSNode* parent; - int child_id; + int child_id; bool operator>(const PqItem& o) const { return lb > o.lb; } }; diff --git a/lib/algos/DumpyOS.hpp b/lib/algos/DumpyOS.hpp index f6f8c60..a4a7886 100644 --- a/lib/algos/DumpyOS.hpp +++ b/lib/algos/DumpyOS.hpp @@ -16,11 +16,11 @@ struct DumpyOSConfig { float fill_upper = 3.0f; // f_high in the paper }; -// Adapted from FADASNode in DumpyOS/include/DataStructures/FADASNode.h +// Adapted from FADASNode in DumpyOS struct DumpyOSNode { std::vector levels; // bits_cardinality[] per segment std::vector sax_word; // SAX word at current bit depth (needed for LB) - std::vector chosen_segs; // chosenSegments; empty ↔ leaf + std::vector chosen_segs; // chosen segments: empty , leaf std::vector children; // 2^|chosen_segs| entries (may be nullptr) std::vector entries; // series indices (leaf only) int n = 0; @@ -33,6 +33,8 @@ class DumpyOS : public SimilaritySearchAlgorithm { using SimilaritySearchAlgorithm::buildIndex; + void setWarpingWindow(int w) { warping_window = w; } + void buildIndex(DataSource* data_source) override; void searchIndex(const float* query, idx_t n_query, idx_t k, idx_t* I, float* D) override; @@ -50,4 +52,4 @@ class DumpyOS : public SimilaritySearchAlgorithm { void destroyTree_(DumpyOSNode* node); }; -} // namespace daisy +} diff --git a/pybinds/setup.cpp b/pybinds/setup.cpp index 634b671..cb56223 100644 --- a/pybinds/setup.cpp +++ b/pybinds/setup.cpp @@ -612,6 +612,7 @@ PYBIND11_MODULE(_core, m) .def(pybind11::init(), "Create a new DumpyOS with the given distance metric") .def(pybind11::init(), "Create a new DumpyOS with the given distance metric and configuration") .def("setNumThreads", &daisy::DumpyOS::setNumThreads, "Set the number of threads") + .def("setWarpingWindow", &daisy::DumpyOS::setWarpingWindow, "Set the warping window size for DTW") .def("buildIndex", [](daisy::DumpyOS &self, pybind11::array_t db) { pybind11::buffer_info buf = db.request(); From dbc4206c3437b1bdf30c89572b6c75c519b69ac2 Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Wed, 1 Jul 2026 15:18:41 +0200 Subject: [PATCH 10/19] Add Fresh algo + tests and demo --- demos/CMakeLists.txt | 31 + demos/demo_Fresh_L2Square.cpp | 42 ++ lib/algos/CMakeLists.txt | 12 + lib/algos/Fresh.cpp | 1307 +++++++++++++++++++++++++++++++++ lib/algos/Fresh.hpp | 175 +++++ lib/daisy.hpp | 1 + lib/isax/iSAXIndex.cpp | 1 + lib/isax/iSAXIndex.hpp | 3 + pybinds/setup.cpp | 71 ++ tests/CMakeLists.txt | 27 + tests/test_Fresh_L2Square.cpp | 45 ++ tests/test_utils.hpp | 11 + 12 files changed, 1726 insertions(+) create mode 100644 demos/demo_Fresh_L2Square.cpp create mode 100644 lib/algos/Fresh.cpp create mode 100644 lib/algos/Fresh.hpp create mode 100644 tests/test_Fresh_L2Square.cpp diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt index f85c069..0d556ad 100644 --- a/demos/CMakeLists.txt +++ b/demos/CMakeLists.txt @@ -496,6 +496,37 @@ if(BUILD_DEMO) message(STATUS "Include directories added for demo_DumpyOS_DTW.") endif() + # ////// FRESH L2Square ////// + if(DEBUG_MSG) + message(STATUS "---") + message(STATUS "## Demo: Fresh L2Square") + message(STATUS "Attempting to add executable: demo_Fresh_L2Square") + endif() + add_executable(demo_Fresh_L2Square demo_Fresh_L2Square.cpp) + if(DEBUG_MSG) + message(STATUS "Executable demo_Fresh_L2Square added.") + endif() + + if(DEBUG_MSG) + message(STATUS "Linking libraries for demo_Fresh_L2Square...") + endif() + target_link_libraries(demo_Fresh_L2Square PRIVATE dino_lib commons_lib) + if(DEBUG_MSG) + message(STATUS "Libraries linked for demo_Fresh_L2Square.") + endif() + + if(DEBUG_MSG) + message(STATUS "Adding include directories for demo_Fresh_L2Square...") + endif() + target_include_directories(demo_Fresh_L2Square PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons + ) + if(DEBUG_MSG) + message(STATUS "Include directories added for demo_Fresh_L2Square.") + endif() + + else() if(DEBUG_MSG) message(STATUS "BUILD_DEMO is FALSE. Demo executables will NOT be built.") diff --git a/demos/demo_Fresh_L2Square.cpp b/demos/demo_Fresh_L2Square.cpp new file mode 100644 index 0000000..7f9c5bf --- /dev/null +++ b/demos/demo_Fresh_L2Square.cpp @@ -0,0 +1,42 @@ +#include "../commons/dataloaders.hpp" +#include "../lib/daisy.hpp" +#include + +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::Fresh fresh_search(daisy::DistanceType::L2_SQUARED); + fresh_search.setNumThreads(4); + + fresh_search.buildIndex(database, n_database, dim); + + daisy::idx_t *I = new daisy::idx_t[n_query * k]; + float *D = new float[n_query * k]; + fresh_search.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/lib/algos/CMakeLists.txt b/lib/algos/CMakeLists.txt index 33bc218..2981ef5 100644 --- a/lib/algos/CMakeLists.txt +++ b/lib/algos/CMakeLists.txt @@ -164,6 +164,18 @@ if(DEBUG_MSG) message(STATUS "DumpyOS.cpp added.") endif() +# ////// FRESH ////// +if(DEBUG_MSG) + message(STATUS "Adding Fresh.cpp to dino_lib sources.") +endif() +target_sources(dino_lib + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/Fresh.cpp +) +if(DEBUG_MSG) + message(STATUS "Fresh.cpp added.") +endif() + # ////// PARIS ////// if(DEBUG_MSG) message(STATUS "Adding ParIS.cpp to dino_lib sources.") diff --git a/lib/algos/Fresh.cpp b/lib/algos/Fresh.cpp new file mode 100644 index 0000000..3e10419 --- /dev/null +++ b/lib/algos/Fresh.cpp @@ -0,0 +1,1307 @@ +#include "Fresh.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include "../isax/iSAXPqueue.hpp" +#include "../isax/iSAXIndex.hpp" +#include "../isax/iSAXSearch.hpp" +#include "../distance_computers/DistanceComputer.hpp" + +namespace daisy +{ + +#define FRESH_TS_GROUP_LENGTH 1000 + +typedef struct { + volatile unsigned long num __attribute__((aligned(64))); +} fresh_next_ts_group_t; + +typedef struct { + volatile unsigned long num __attribute__((aligned(64))); +} fresh_next_ts_in_group_t; + +typedef struct FRESH_buffer_data { + isax_index *index; + int workernumber; + int total_workernumber; + unsigned long ts_num; + unsigned long *shared_start_number; + pthread_barrier_t *wait_summaries_to_compute; + int *node_counter; + float *rawfile; + volatile unsigned char *block_processed; + volatile unsigned char **group_processed; + volatile unsigned char *ts_processed; + fresh_next_ts_group_t *next_ts_group_read_in_block; + fresh_next_ts_in_group_t **next_ts_read_in_group; + fresh_next_ts_in_group_t **next_ts_read_in_group_fai; + volatile unsigned char *all_blocks_processed; + volatile unsigned char *all_RecBufs_processed; + volatile unsigned char *block_helper_exist; + volatile unsigned char **group_helpers_exist; + unsigned long total_blocks; + unsigned long total_groups; + unsigned long read_block_length; +} FRESH_buffer_data; + +static volatile int g_fresh_query_id = 0; + +// defined in Messi.cpp, same namespace +void calculate_node2_topk_inmemory(isax_index *index, isax_node *node, ts_type *query, ts_type *paa, + pqueue_bsf *pq_bsf, pthread_rwlock_t *lock_queue, float *rawfile); +void calculate_node_DTW2knn_inmemory(isax_index *index, isax_node *node, ts_type *query, + float *uo, float *lo, ts_type *paa, ts_type *paaU, ts_type *paaL, + float bsf, int warpWind, pqueue_bsf *pq_bsf, + pthread_rwlock_t *lock_queue, float *rawfile); + +static parallel_first_buffer_layer_lf *initialize_pRecBuf_lf(int initial_buffer_size, int number_of_buffers, + int max_total_size, isax_index *index) +{ + parallel_first_buffer_layer_lf *fbl = (parallel_first_buffer_layer_lf *)malloc(sizeof(parallel_first_buffer_layer_lf)); + fbl->number_of_buffers = number_of_buffers; + fbl->initial_buffer_size = initial_buffer_size; + fbl->max_total_size = max_total_size; + fbl->current_record_index = 0; + fbl->soft_buffers = (parallel_fbl_soft_buffer_lf *)calloc(number_of_buffers, sizeof(parallel_fbl_soft_buffer_lf)); + return fbl; +} + +static isax_node *insert_to_pRecBuf_lock_free(parallel_first_buffer_layer_lf *fbl, sax_type *sax, + file_position_type *pos, root_mask_type mask, + isax_index *index, int workernumber, int total_workernumber) +{ + parallel_fbl_soft_buffer_lf *current_buffer = &fbl->soft_buffers[(int)mask]; + current_buffer->mask = mask; + + if (!current_buffer->initialized) + { + int *tmp_max_buffer_size = (int *)calloc(total_workernumber, sizeof(int)); + int *tmp_buffer_size = (int *)calloc(total_workernumber, sizeof(int)); + sax_type **tmp_sax_records = (sax_type **)calloc(total_workernumber, sizeof(sax_type *)); + file_position_type **tmp_pos_records = (file_position_type **)calloc(total_workernumber, sizeof(file_position_type *)); + + if (!current_buffer->max_buffer_size && !FRESH_CASPTR(¤t_buffer->max_buffer_size, NULL, tmp_max_buffer_size)) + free(tmp_max_buffer_size); + if (!current_buffer->buffer_size && !FRESH_CASPTR(¤t_buffer->buffer_size, NULL, tmp_buffer_size)) + free(tmp_buffer_size); + if (!current_buffer->sax_records && !FRESH_CASPTR(¤t_buffer->sax_records, NULL, tmp_sax_records)) + free(tmp_sax_records); + if (!current_buffer->pos_records && !FRESH_CASPTR(¤t_buffer->pos_records, NULL, tmp_pos_records)) + free(tmp_pos_records); + + current_buffer->mask = mask; + if (!current_buffer->initialized) + current_buffer->initialized = 1; + } + + if (current_buffer->buffer_size[workernumber] >= current_buffer->max_buffer_size[workernumber]) + { + if (current_buffer->max_buffer_size[workernumber] == 0) + { + current_buffer->max_buffer_size[workernumber] = fbl->initial_buffer_size; + current_buffer->sax_records[workernumber] = (sax_type *)malloc(index->settings->sax_byte_size * + current_buffer->max_buffer_size[workernumber]); + current_buffer->pos_records[workernumber] = (file_position_type *)malloc(index->settings->position_byte_size * + current_buffer->max_buffer_size[workernumber]); + } + else + { + current_buffer->max_buffer_size[workernumber] *= BUFFER_REALLOCATION_RATE; + current_buffer->sax_records[workernumber] = (sax_type *)realloc(current_buffer->sax_records[workernumber], + index->settings->sax_byte_size * current_buffer->max_buffer_size[workernumber]); + current_buffer->pos_records[workernumber] = (file_position_type *)realloc(current_buffer->pos_records[workernumber], + index->settings->position_byte_size * current_buffer->max_buffer_size[workernumber]); + } + } + + if (current_buffer->sax_records[workernumber] == NULL || current_buffer->pos_records[workernumber] == NULL) + { + fprintf(stderr, "error: Could not allocate memory in FBL."); + return NULL; + } + + int current_buffer_number = current_buffer->buffer_size[workernumber]; + file_position_type *filepointer = (file_position_type *)current_buffer->pos_records[workernumber]; + sax_type *saxpointer = (sax_type *)current_buffer->sax_records[workernumber]; + + memcpy(&saxpointer[current_buffer_number * index->settings->paa_segments], sax, index->settings->sax_byte_size); + memcpy(&filepointer[current_buffer_number], pos, index->settings->position_byte_size); + + (current_buffer->buffer_size[workernumber])++; + return (isax_node *)current_buffer->node; +} + +static void store_isax_in_pRecBuf_fresh(FRESH_buffer_data *input_data, isax_index *index, + unsigned long ts_id, sax_type *sax) +{ + file_position_type pos; + + if (sax_from_ts((ts_type *)&input_data->rawfile[ts_id * index->settings->timeseries_size], + sax, + index->settings->ts_values_per_paa_segment, + index->settings->paa_segments, + index->settings->sax_alphabet_cardinality, + index->settings->sax_bit_cardinality) == SUCCESS) + { + pos = (file_position_type)(ts_id * index->settings->timeseries_size); + + root_mask_type first_bit_mask = 0x00; + CREATE_MASK(first_bit_mask, index, sax); + + insert_to_pRecBuf_lock_free((parallel_first_buffer_layer_lf *)(index->fbl), + sax, &pos, first_bit_mask, index, + input_data->workernumber, input_data->total_workernumber); + } + else + { + fprintf(stderr, "error: cannot insert record in index, since sax representation failed to be created"); + } +} + +static void scan_for_unprocessed_ts_fresh(FRESH_buffer_data *input_data, isax_index *index, + unsigned long ts_start, unsigned long ts_end, + volatile unsigned char *stop, sax_type *sax) +{ + for (unsigned long ts_id = ts_start; ts_id < ts_end && !(*stop); ts_id++) + { + if (!input_data->ts_processed[ts_id]) + { + store_isax_in_pRecBuf_fresh(input_data, index, ts_id, sax); + if (!input_data->ts_processed[ts_id]) + input_data->ts_processed[ts_id] = 1; + } + } +} + +static void process_group_fresh(unsigned long ts_group, unsigned long block_num, + unsigned long ts_group_start, unsigned long ts_group_end, + FRESH_buffer_data *input_data, isax_index *index, + char is_helper, sax_type *sax) +{ + unsigned long ts_id, tmp; + + while (!input_data->group_processed[block_num][ts_group]) + { + if (!input_data->group_helpers_exist[block_num][ts_group]) + { + tmp = input_data->next_ts_read_in_group[block_num][ts_group].num; + input_data->next_ts_read_in_group[block_num][ts_group].num = tmp + 1; + } + else + { + if (!input_data->next_ts_read_in_group_fai[block_num][ts_group].num && + input_data->next_ts_read_in_group[block_num][ts_group].num) + { + FRESH_CASULONG(&input_data->next_ts_read_in_group_fai[block_num][ts_group].num, + 0, input_data->next_ts_read_in_group[block_num][ts_group].num); + } + tmp = __sync_fetch_and_add(&input_data->next_ts_read_in_group_fai[block_num][ts_group].num, 1); + } + + ts_id = ts_group_start + tmp; + if (ts_id >= ts_group_end) + break; + + if (!input_data->ts_processed[ts_id]) + { + store_isax_in_pRecBuf_fresh(input_data, index, ts_id, sax); + if (!input_data->ts_processed[ts_id]) + input_data->ts_processed[ts_id] = 1; + } + } + + if ((is_helper || input_data->group_helpers_exist[block_num][ts_group]) && + !input_data->group_processed[block_num][ts_group]) + { + scan_for_unprocessed_ts_fresh(input_data, index, ts_group_start, ts_group_end, + &input_data->group_processed[block_num][ts_group], sax); + } + + if (!input_data->group_processed[block_num][ts_group]) + input_data->group_processed[block_num][ts_group] = 1; +} + +static void scan_for_unprocessed_groups_fresh(unsigned long block_num, FRESH_buffer_data *input_data, + isax_index *index, unsigned long total_groups_in_block, + unsigned long my_ts_start, unsigned long my_ts_end, + volatile unsigned char *stop, sax_type *sax) +{ + for (unsigned long ts_group = 0; ts_group < total_groups_in_block && !*stop; ts_group++) + { + unsigned long ts_group_start = my_ts_start + ts_group * FRESH_TS_GROUP_LENGTH; + unsigned long ts_group_end = (ts_group == total_groups_in_block - 1) + ? my_ts_end + : ts_group_start + FRESH_TS_GROUP_LENGTH; + + if (!input_data->group_helpers_exist[block_num][ts_group]) + input_data->group_helpers_exist[block_num][ts_group] = 1; + + process_group_fresh(ts_group, block_num, ts_group_start, ts_group_end, input_data, index, 1, sax); + } + + if (!input_data->block_processed[block_num]) + input_data->block_processed[block_num] = 1; +} + +static void process_block_fresh(unsigned long block_num, unsigned long total_blocks, + unsigned long total_ts_num, FRESH_buffer_data *input_data, + isax_index *index, char is_helper, sax_type *sax) +{ + unsigned long my_ts_start = block_num * input_data->read_block_length; + unsigned long my_ts_end = (block_num == total_blocks - 1) + ? total_ts_num + : my_ts_start + input_data->read_block_length; + + unsigned long total_groups_in_block = (my_ts_end - my_ts_start) / FRESH_TS_GROUP_LENGTH; + if (total_groups_in_block * FRESH_TS_GROUP_LENGTH < my_ts_end - my_ts_start) + total_groups_in_block++; + + unsigned long prev_group_id = 0; + char helpers_exist = 0; + + while (!input_data->block_processed[block_num]) + { + unsigned long ts_group; + + if (!input_data->block_helper_exist[block_num]) + { + ts_group = input_data->next_ts_group_read_in_block[block_num].num; + input_data->next_ts_group_read_in_block[block_num].num = ts_group + 1; + } + else + { + ts_group = __sync_fetch_and_add(&input_data->next_ts_group_read_in_block[block_num].num, 1); + } + + if (ts_group > prev_group_id + 1) + helpers_exist = 1; + + if (ts_group >= total_groups_in_block) + break; + + unsigned long ts_group_start = my_ts_start + ts_group * FRESH_TS_GROUP_LENGTH; + unsigned long ts_group_end = (ts_group == total_groups_in_block - 1) + ? my_ts_end + : ts_group_start + FRESH_TS_GROUP_LENGTH; + + if (is_helper && !input_data->group_helpers_exist[block_num][ts_group]) + input_data->group_helpers_exist[block_num][ts_group] = 1; + + process_group_fresh(ts_group, block_num, ts_group_start, ts_group_end, input_data, index, is_helper, sax); + + prev_group_id = ts_group; + } + + if ((is_helper || helpers_exist) && !input_data->block_processed[block_num]) + scan_for_unprocessed_groups_fresh(block_num, input_data, index, total_groups_in_block, + my_ts_start, my_ts_end, &input_data->block_processed[block_num], sax); + + if (!input_data->block_processed[block_num]) + input_data->block_processed[block_num] = 1; +} + +static void scan_for_unprocessed_blocks_fresh(FRESH_buffer_data *input_data, isax_index *index, + unsigned long total_blocks, sax_type *sax) +{ + unsigned long total_ts_num = input_data->ts_num; + unsigned long my_id = input_data->workernumber; + unsigned long start_block_num = (my_id + 1) % total_blocks; + unsigned long block_num = start_block_num; + + do + { + if (*input_data->all_blocks_processed) + break; + + if (input_data->block_processed[block_num]) + { + block_num = (block_num + 1) % total_blocks; + continue; + } + + if (!input_data->block_helper_exist[block_num]) + input_data->block_helper_exist[block_num] = 1; + + process_block_fresh(block_num, total_blocks, total_ts_num, input_data, index, 1, sax); + + block_num = (block_num + 1) % total_blocks; + } while (block_num != start_block_num); + + if (!*input_data->all_blocks_processed) + *input_data->all_blocks_processed = 1; +} + +static void *indexCreationWorkerFresh(void *transferdata) +{ + FRESH_buffer_data *input_data = (FRESH_buffer_data *)transferdata; + isax_index *index = input_data->index; + + unsigned long total_ts_num = input_data->ts_num; + unsigned long total_blocks = input_data->total_blocks; + + sax_type *sax = (sax_type *)malloc(sizeof(sax_type) * index->settings->paa_segments); + + // Phase 1: summarization — FAI on blocks, Refresh within each block + unsigned long block_num; + while (!*input_data->all_blocks_processed) + { + block_num = __sync_fetch_and_add(input_data->shared_start_number, 1); + if (block_num >= total_blocks) + break; + + process_block_fresh(block_num, total_blocks, total_ts_num, input_data, index, 0, sax); + } + + scan_for_unprocessed_blocks_fresh(input_data, index, total_blocks, sax); + + pthread_barrier_wait(input_data->wait_summaries_to_compute); + + // Phase 2: tree population — FAI on FBL slots, serial add_record_to_node_inmemory + isax_node_record *r = (isax_node_record *)malloc(sizeof(isax_node_record)); + + while (1) + { + int j = __sync_fetch_and_add(input_data->node_counter, 1); + if (j >= index->fbl->number_of_buffers) + break; + + parallel_fbl_soft_buffer_lf *current_fbl_node = + &((parallel_first_buffer_layer_lf *)(index->fbl))->soft_buffers[j]; + + if (!current_fbl_node->initialized) + continue; + + bool have_record = false; + for (int k = 0; k < input_data->total_workernumber; k++) + { + if (current_fbl_node->buffer_size[k] > 0) + have_record = true; + for (int i = 0; i < current_fbl_node->buffer_size[k]; i++) + { + r->sax = (sax_type *)&((current_fbl_node->sax_records[k])[i * index->settings->paa_segments]); + r->position = (file_position_type *)&((file_position_type *)(current_fbl_node->pos_records[k]))[i]; + r->insertion_mode = (insertion_mode)(NO_TMP | PARTIAL); + + if (!current_fbl_node->node) + { + isax_node *root = isax_root_node_init(current_fbl_node->mask, index->settings->initial_leaf_buffer_size); + if (!FRESH_CASPTR(¤t_fbl_node->node, NULL, root)) + free(root); + } + + add_record_to_node_inmemory(index, current_fbl_node->node, r, 1); + } + } + + if (have_record) + flush_subtree_leaf_buffers_inmemory(index, current_fbl_node->node); + } + + free(r); + free(sax); + return NULL; +} + +// ── Phase 3: lock-free candidate queue (array list + sorted array) ──────────── + +static void init_array_list_fresh(fresh_array_list_t *list, int size) +{ + list->Top = (fresh_array_list_node_t *)malloc(sizeof(fresh_array_list_node_t)); + list->Top->data = (fresh_array_element_t *)calloc(size, sizeof(fresh_array_element_t)); + list->Top->num_node = 0; + list->Top->next = NULL; + list->element_size = size; +} + +static fresh_array_list_node_t *add_array_fresh(fresh_array_list_node_t *last, fresh_array_list_t *list) +{ + fresh_array_list_node_t *tmp = (fresh_array_list_node_t *)malloc(sizeof(fresh_array_list_node_t)); + tmp->data = (fresh_array_element_t *)calloc(list->element_size, sizeof(fresh_array_element_t)); + tmp->num_node = last->num_node + 1; + tmp->next = last; + + if (!FRESH_CASPTR(&list->Top, last, tmp)) + { + free(tmp->data); + free(tmp); + tmp = list->Top; + while (tmp->num_node > last->num_node + 1) + tmp = tmp->next; + } + return tmp; +} + +static fresh_array_element_t *get_element_at_fresh(unsigned long position, fresh_array_list_t *list) +{ + fresh_array_list_node_t *array_node; + int array_num_node = position / list->element_size; + fresh_array_list_node_t *last = list->Top; + + if (array_num_node == last->num_node) + array_node = last; + else if (array_num_node > last->num_node) + array_node = add_array_fresh(last, list); + else + { + array_node = list->Top; + while (array_node->num_node > array_num_node) + array_node = array_node->next; + } + return &array_node->data[position % list->element_size]; +} + +static void add_to_array_data_lf(float *paa, isax_node *node, isax_index *index, float bsf, + fresh_array_list_t *array_lists, int *tnumber, + volatile unsigned long *next_queue_data_pos, int n_queues, + const int query_id) +{ + if (node->processed == (unsigned long)query_id) + return; + + float distance = minidist_paa_to_isax(paa, node->isax_values, + node->isax_cardinalities, + index->settings->sax_bit_cardinality, + index->settings->sax_alphabet_cardinality, + index->settings->paa_segments, + MINVAL, MAXVAL, + index->settings->mindist_sqrt); + + if (distance < bsf && node->processed < (unsigned long)query_id) + { + if (node->is_leaf) + { + unsigned long queue_data_pos = __sync_fetch_and_add(&next_queue_data_pos[*tnumber], 1); + fresh_array_element_t *array_elem = get_element_at_fresh(queue_data_pos, &array_lists[*tnumber]); + array_elem->distance = distance; + array_elem->node = node; + *tnumber = (*tnumber + 1) % n_queues; + } + else + { + if (node->left_child->isax_cardinalities != NULL && node->left_child->processed < (unsigned long)query_id) + add_to_array_data_lf(paa, node->left_child, index, bsf, array_lists, tnumber, next_queue_data_pos, n_queues, query_id); + if (node->right_child->isax_cardinalities != NULL && node->right_child->processed < (unsigned long)query_id) + add_to_array_data_lf(paa, node->right_child, index, bsf, array_lists, tnumber, next_queue_data_pos, n_queues, query_id); + } + } + + if (node->processed < (unsigned long)query_id) + node->processed = query_id; +} + +static int compare_sorted_array_items(const void *a, const void *b) +{ + float dist_a = ((fresh_array_element_t *)a)->distance; + float dist_b = ((fresh_array_element_t *)b)->distance; + if (dist_a < dist_b) return -1; + if (dist_a > dist_b) return 1; + return 0; +} + +static void create_sorted_array_from_data_queue(int pq_id, fresh_array_list_t *array_lists, + volatile unsigned long *next_queue_data_pos, + fresh_sorted_array_t *volatile *sorted_arrays) +{ + fresh_sorted_array_t *local_sa = (fresh_sorted_array_t *)malloc(sizeof(fresh_sorted_array_t)); + unsigned long array_elements = next_queue_data_pos[pq_id]; + local_sa->data = (fresh_array_element_t *)malloc(array_elements * sizeof(fresh_array_element_t)); + + size_t j = 0; + int array_buckets_num = array_lists[pq_id].Top->num_node + 1; + unsigned long array_elements_traversed = 0; + fresh_array_list_node_t *bucket = array_lists[pq_id].Top; + for (int i = 0; i < array_buckets_num; i++, bucket = bucket->next) + { + for (int k = 0; k < array_lists[pq_id].element_size && + array_elements_traversed < array_elements && !sorted_arrays[pq_id]; + k++, array_elements_traversed++) + { + if (bucket->data[k].node && bucket->data[k].distance > 0) + { + local_sa->data[j].node = bucket->data[k].node; + local_sa->data[j].distance = bucket->data[k].distance; + j++; + } + } + } + + local_sa->num_elements = j; + if (!sorted_arrays[pq_id]) + qsort(local_sa->data, j, sizeof(fresh_array_element_t), compare_sorted_array_items); + + if (sorted_arrays[pq_id] || !FRESH_CASPTR(&sorted_arrays[pq_id], NULL, local_sa)) + { + free(local_sa->data); + free(local_sa); + } +} + +// ── Phase 4: refinement (help sorted arrays) ────────────────────────────────── + +static int process_sorted_array_element_L2(fresh_array_element_t *n, FRESH_workerdata *wd) +{ + float bsfdistance = wd->pq_bsf->knn[wd->pq_bsf->k - 1]; + if (n->distance > bsfdistance || n->distance > wd->minimum_distance) + return 0; + if (n->node->is_leaf && n->distance >= 0) + { + calculate_node2_topk_inmemory(wd->index, n->node, wd->ts, wd->paa, wd->pq_bsf, wd->lock_bsf, wd->rawfile); + n->distance = -1; + } + return 1; +} + +static int process_sorted_array_element_DTW(fresh_array_element_t *n, FRESH_workerdata *wd) +{ + float bsfdistance = wd->pq_bsf->knn[wd->pq_bsf->k - 1]; + if (n->distance > bsfdistance || n->distance > wd->minimum_distance) + return 0; + if (n->node->is_leaf && n->distance >= 0) + { + calculate_node_DTW2knn_inmemory(wd->index, n->node, wd->ts, wd->uo, wd->lo, + wd->paa, wd->paaU, wd->paaL, bsfdistance, + wd->warpWind, wd->pq_bsf, wd->lock_bsf, wd->rawfile); + n->distance = -1; + } + return 1; +} + +static void help_sorted_array_L2(FRESH_workerdata *wd, fresh_sorted_array_t *sa, int pq_id) +{ + size_t pq_size = sa->num_elements; + if (pq_size != 0) + { + int element_id; + while ((element_id = (int)__sync_fetch_and_add(&wd->sorted_array_FAI_counter[pq_id], 1)) < (int)pq_size) + { + fresh_array_element_t *n = &sa->data[element_id]; + if (!process_sorted_array_element_L2(n, wd)) + break; + } + for (int i = 0; i < (int)pq_size && !wd->queue_finished[pq_id]; i++) + { + fresh_array_element_t *n = &sa->data[i]; + if (n->distance >= 0 && !process_sorted_array_element_L2(n, wd)) + break; + } + } + if (!wd->queue_finished[pq_id]) + wd->queue_finished[pq_id] = 1; +} + +static void help_sorted_array_DTW(FRESH_workerdata *wd, fresh_sorted_array_t *sa, int pq_id) +{ + size_t pq_size = sa->num_elements; + if (pq_size != 0) + { + int element_id; + while ((element_id = (int)__sync_fetch_and_add(&wd->sorted_array_FAI_counter[pq_id], 1)) < (int)pq_size) + { + fresh_array_element_t *n = &sa->data[element_id]; + if (!process_sorted_array_element_DTW(n, wd)) + break; + } + for (int i = 0; i < (int)pq_size && !wd->queue_finished[pq_id]; i++) + { + fresh_array_element_t *n = &sa->data[i]; + if (n->distance >= 0 && !process_sorted_array_element_DTW(n, wd)) + break; + } + } + if (!wd->queue_finished[pq_id]) + wd->queue_finished[pq_id] = 1; +} + +// ── Search workers ───────────────────────────────────────────────────────────── + +void *FRESH_topk_search_worker_L2Squared(void *rfdata) +{ + FRESH_workerdata *wd = (FRESH_workerdata *)rfdata; + isax_index *index = wd->index; + ts_type *paa = wd->paa; + float bsfdistance = wd->pq_bsf->knn[wd->pq_bsf->k - 1]; + const int query_id = wd->query_id; + int tnumber = rand() % wd->n_queues; + + // A. populate arrays + while (1) + { + int current_root_node_number = __sync_fetch_and_add(wd->node_counter, 1); + if (current_root_node_number >= wd->amountnode) + break; + add_to_array_data_lf(paa, wd->nodelist[current_root_node_number], index, bsfdistance, + wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); + } + + // A.1. help unprocessed subtrees + for (int i = 0; i < wd->amountnode; i++) + { + isax_node *node = wd->nodelist[i]; + if (node->processed == (unsigned long)query_id) + continue; + add_to_array_data_lf(paa, node, index, bsfdistance, + wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); + } + + // A.2. create my sorted array + int my_pq = wd->workernumber % wd->n_queues; + if (!wd->sorted_arrays[my_pq] && wd->next_queue_data_pos[my_pq]) + create_sorted_array_from_data_queue(my_pq, wd->array_lists, wd->next_queue_data_pos, wd->sorted_arrays); + + // B. process my sorted array + fresh_sorted_array_t *my_array = wd->sorted_arrays[my_pq]; + if (my_array && !wd->queue_finished[my_pq]) + help_sorted_array_L2(wd, my_array, my_pq); + + // A.3. help create other sorted arrays + for (int i = 0; i < wd->n_queues; i++) + { + if (!wd->sorted_arrays[i] && wd->next_queue_data_pos[i]) + create_sorted_array_from_data_queue(i, wd->array_lists, wd->next_queue_data_pos, wd->sorted_arrays); + } + + // B.3. help process other sorted arrays + for (int i = 0; i < wd->n_queues; i++) + { + if (wd->queue_finished[i] || !wd->sorted_arrays[i]) + continue; + if (!wd->helper_queue_exist[i]) + wd->helper_queue_exist[i] = 1; + help_sorted_array_L2(wd, wd->sorted_arrays[i], i); + } + + return nullptr; +} + +void *FRESH_topk_search_worker_DTW(void *rfdata) +{ + FRESH_workerdata *wd = (FRESH_workerdata *)rfdata; + isax_index *index = wd->index; + ts_type *paa = wd->paa; + float bsfdistance = wd->pq_bsf->knn[wd->pq_bsf->k - 1]; + const int query_id = wd->query_id; + int tnumber = rand() % wd->n_queues; + + // A. populate arrays + while (1) + { + int current_root_node_number = __sync_fetch_and_add(wd->node_counter, 1); + if (current_root_node_number >= wd->amountnode) + break; + add_to_array_data_lf(paa, wd->nodelist[current_root_node_number], index, bsfdistance, + wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); + } + + // A.1. help unprocessed subtrees + for (int i = 0; i < wd->amountnode; i++) + { + isax_node *node = wd->nodelist[i]; + if (node->processed == (unsigned long)query_id) + continue; + add_to_array_data_lf(paa, node, index, bsfdistance, + wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); + } + + // A.2. create my sorted array + int my_pq = wd->workernumber % wd->n_queues; + if (!wd->sorted_arrays[my_pq] && wd->next_queue_data_pos[my_pq]) + create_sorted_array_from_data_queue(my_pq, wd->array_lists, wd->next_queue_data_pos, wd->sorted_arrays); + + // B. process my sorted array + fresh_sorted_array_t *my_array = wd->sorted_arrays[my_pq]; + if (my_array && !wd->queue_finished[my_pq]) + help_sorted_array_DTW(wd, my_array, my_pq); + + // A.3. help create other sorted arrays + for (int i = 0; i < wd->n_queues; i++) + { + if (!wd->sorted_arrays[i] && wd->next_queue_data_pos[i]) + create_sorted_array_from_data_queue(i, wd->array_lists, wd->next_queue_data_pos, wd->sorted_arrays); + } + + // B.3. help process other sorted arrays + for (int i = 0; i < wd->n_queues; i++) + { + if (wd->queue_finished[i] || !wd->sorted_arrays[i]) + continue; + if (!wd->helper_queue_exist[i]) + wd->helper_queue_exist[i] = 1; + help_sorted_array_DTW(wd, wd->sorted_arrays[i], i); + } + + return nullptr; +} + +// ── Fresh class ──────────────────────────────────────────────────────────────── + +Fresh::Fresh(DistanceType distance_type) + : Fresh(distance_type, FreshConfig{}) +{ +} + +Fresh::Fresh(DistanceType distance_type, const FreshConfig &config) + : SimilaritySearchAlgorithm(distance_type) +{ + this->search_workers = config.search_workers; + this->index_workers = config.index_workers; + this->warping_window = config.warping_window; + this->leaf_size = config.leaf_size; + this->paa_segments = config.paa_segments; +} + +void Fresh::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(); + } + + data_source->reset(); + const float *raw = data_source->rawPointer(); + if (raw != nullptr) + { + this->database = const_cast(raw); + this->owns_database = false; + } + else + { + this->database = new float[this->n_database * this->dim]; + float *record = new float[this->dim]; + idx_t idx = 0; + while (data_source->nextRecord(record)) + { + std::copy(record, record + this->dim, this->database + idx * this->dim); + idx++; + } + delete[] record; + this->owns_database = true; + } + + this->index_settings = isax_index_settings_init("", + this->dim, + this->paa_segments, + this->sax_cardinality, + 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->index = isax_index_init_inmemory(this->index_settings); + isax_index *index = this->index; + index->sax_file = NULL; + + unsigned long shared_start_number = 0; + int node_counter = 0; + + int n_buffers = (int)pow(2, index->settings->paa_segments); + destroy_fbl(index->fbl); + index->fbl = (first_buffer_layer *)initialize_pRecBuf_lf( + index->settings->initial_fbl_buffer_size, n_buffers, + index->settings->max_total_buffer_size + DISK_BUFFER_SIZE * (PROGRESS_CALCULATE_THREAD_NUMBER - 1), + index); + + unsigned long total_blocks = this->n_database / this->read_block_length; + if (this->read_block_length * total_blocks < this->n_database) + total_blocks++; + + unsigned long total_groups = this->read_block_length / FRESH_TS_GROUP_LENGTH; + if (FRESH_TS_GROUP_LENGTH * total_groups < (unsigned long)this->read_block_length) + total_groups++; + + volatile unsigned char *block_processed = (volatile unsigned char *)calloc(total_blocks, sizeof(unsigned char)); + volatile unsigned char **group_processed = (volatile unsigned char **)malloc(total_blocks * sizeof(unsigned char *)); + fresh_next_ts_group_t *next_ts_group_read_in_block = (fresh_next_ts_group_t *)calloc(total_blocks, sizeof(fresh_next_ts_group_t)); + fresh_next_ts_in_group_t **next_ts_read_in_group = (fresh_next_ts_in_group_t **)malloc(total_blocks * sizeof(fresh_next_ts_in_group_t *)); + fresh_next_ts_in_group_t **next_ts_read_in_group_fai = (fresh_next_ts_in_group_t **)malloc(total_blocks * sizeof(fresh_next_ts_in_group_t *)); + volatile unsigned char *block_helper_exist = (volatile unsigned char *)calloc(total_blocks, sizeof(unsigned char)); + volatile unsigned char **group_helpers_exist = (volatile unsigned char **)malloc(total_blocks * sizeof(unsigned char *)); + + for (unsigned long i = 0; i < total_blocks; i++) + { + group_processed[i] = (volatile unsigned char *)calloc(total_groups, sizeof(unsigned char)); + next_ts_read_in_group[i] = (fresh_next_ts_in_group_t *)calloc(total_groups, sizeof(fresh_next_ts_in_group_t)); + next_ts_read_in_group_fai[i] = (fresh_next_ts_in_group_t *)calloc(total_groups, sizeof(fresh_next_ts_in_group_t)); + group_helpers_exist[i] = (volatile unsigned char *)calloc(total_groups, sizeof(unsigned char)); + } + + volatile unsigned char *ts_processed = (volatile unsigned char *)calloc(this->n_database, sizeof(unsigned char)); + volatile unsigned char all_blocks_processed = 0; + volatile unsigned char all_RecBufs_processed = 0; + + index->sax_cache = (sax_type *)malloc(sizeof(sax_type) * index->settings->paa_segments * this->n_database); + + pthread_barrier_t wait_summaries_to_compute; + pthread_barrier_init(&wait_summaries_to_compute, NULL, this->index_workers); + + pthread_t threadid[this->index_workers]; + FRESH_buffer_data *input_data = (FRESH_buffer_data *)malloc(sizeof(FRESH_buffer_data) * this->index_workers); + + for (int i = 0; i < this->index_workers; i++) + { + input_data[i].index = index; + input_data[i].workernumber = i; + input_data[i].total_workernumber = this->index_workers; + input_data[i].ts_num = this->n_database; + input_data[i].shared_start_number = &shared_start_number; + input_data[i].wait_summaries_to_compute = &wait_summaries_to_compute; + input_data[i].node_counter = &node_counter; + input_data[i].rawfile = this->database; + input_data[i].block_processed = block_processed; + input_data[i].group_processed = group_processed; + input_data[i].ts_processed = ts_processed; + input_data[i].next_ts_group_read_in_block = next_ts_group_read_in_block; + input_data[i].next_ts_read_in_group = next_ts_read_in_group; + input_data[i].next_ts_read_in_group_fai = next_ts_read_in_group_fai; + input_data[i].all_blocks_processed = &all_blocks_processed; + input_data[i].all_RecBufs_processed = &all_RecBufs_processed; + input_data[i].block_helper_exist = block_helper_exist; + input_data[i].group_helpers_exist = group_helpers_exist; + input_data[i].total_blocks = total_blocks; + input_data[i].total_groups = total_groups; + input_data[i].read_block_length = this->read_block_length; + } + + for (int i = 0; i < this->index_workers; i++) + pthread_create(&threadid[i], NULL, indexCreationWorkerFresh, (void *)&input_data[i]); + for (int i = 0; i < this->index_workers; i++) + pthread_join(threadid[i], NULL); + + __sync_fetch_and_add(&(index->total_records), this->n_database); + index->sax_cache_size = index->total_records; + fprintf(stderr, ">>> Finished indexing\n"); + + pthread_barrier_destroy(&wait_summaries_to_compute); + + for (unsigned long i = 0; i < total_blocks; i++) + { + free((void *)group_processed[i]); + free(next_ts_read_in_group[i]); + free(next_ts_read_in_group_fai[i]); + free((void *)group_helpers_exist[i]); + } + free((void *)block_processed); + free(group_processed); + free(next_ts_group_read_in_block); + free(next_ts_read_in_group); + free(next_ts_read_in_group_fai); + free((void *)block_helper_exist); + free(group_helpers_exist); + free((void *)ts_processed); + free(input_data); +} + +pqueue_bsf Fresh::FRESH_search_topk_L2Squared(ts_type *ts, ts_type *paa, node_list *nodelist, idx_t k) +{ + pqueue_bsf *pq_bsf = pqueue_bsf_init(k); + + approximate_topk_inmemory(ts, paa, index, pq_bsf, this->database); + this->minimum_distance = pq_bsf->knn[k - 1]; + + if (this->minimum_distance == FLT_MAX || min_checked_leaves > 1) + { + refine_topk_answer_inmemory(ts, paa, index, pq_bsf, this->minimum_distance, this->min_checked_leaves, this->database); + this->minimum_distance = pq_bsf->knn[k - 1]; + } + + int n_queues = (this->search_workers == 1) ? 1 : this->search_workers / 2; + int query_id = __sync_add_and_fetch(&g_fresh_query_id, 1); + + fresh_array_list_t *array_lists = (fresh_array_list_t *)malloc(sizeof(fresh_array_list_t) * n_queues); + for (int i = 0; i < n_queues; i++) + init_array_list_fresh(&array_lists[i], index->settings->root_nodes_size / n_queues); + + fresh_sorted_array_t *volatile *sorted_arrays = (fresh_sorted_array_t *volatile *)calloc(n_queues, sizeof(fresh_sorted_array_t *)); + volatile unsigned long *next_queue_data_pos = (volatile unsigned long *)calloc(n_queues, sizeof(unsigned long)); + volatile unsigned char *queue_finished = (volatile unsigned char *)calloc(n_queues, sizeof(unsigned char)); + volatile unsigned char *helper_queue_exist = (volatile unsigned char *)calloc(n_queues, sizeof(unsigned char)); + volatile unsigned long *sorted_array_FAI_counter = (volatile unsigned long *)calloc(n_queues, sizeof(unsigned long)); + + volatile int node_counter = 0; + + pthread_rwlock_t lock_bsf = PTHREAD_RWLOCK_INITIALIZER; + pthread_t threadid[this->search_workers]; + FRESH_workerdata workerdata[this->search_workers]; + + for (int i = 0; i < this->search_workers; i++) + { + workerdata[i].paa = paa; + workerdata[i].ts = ts; + workerdata[i].index = index; + workerdata[i].minimum_distance = this->minimum_distance; + workerdata[i].pq_bsf = pq_bsf; + workerdata[i].lock_bsf = &lock_bsf; + workerdata[i].nodelist = nodelist->nlist; + workerdata[i].amountnode = nodelist->node_amount; + workerdata[i].node_counter = &node_counter; + workerdata[i].array_lists = array_lists; + workerdata[i].sorted_arrays = sorted_arrays; + workerdata[i].next_queue_data_pos = next_queue_data_pos; + workerdata[i].queue_finished = queue_finished; + workerdata[i].helper_queue_exist = helper_queue_exist; + workerdata[i].sorted_array_FAI_counter = sorted_array_FAI_counter; + workerdata[i].sorted_array_counter = nullptr; + workerdata[i].fai_queue_counters = nullptr; + workerdata[i].queue_bsf_distance = nullptr; + workerdata[i].workernumber = i; + workerdata[i].n_queues = n_queues; + workerdata[i].query_id = query_id; + workerdata[i].rawfile = this->database; + } + + for (int i = 0; i < this->search_workers; i++) + pthread_create(&threadid[i], NULL, FRESH_topk_search_worker_L2Squared, (void *)&workerdata[i]); + for (int i = 0; i < this->search_workers; i++) + pthread_join(threadid[i], NULL); + + this->minimum_distance = pq_bsf->knn[k - 1]; + + for (int i = 0; i < n_queues; i++) + { + if (sorted_arrays[i]) + { + free(sorted_arrays[i]->data); + free(sorted_arrays[i]); + } + fresh_array_list_node_t *node = array_lists[i].Top; + while (node) + { + fresh_array_list_node_t *next = node->next; + free(node->data); + free(node); + node = next; + } + } + free(array_lists); + free((void *)sorted_arrays); + free((void *)next_queue_data_pos); + free((void *)queue_finished); + free((void *)helper_queue_exist); + free((void *)sorted_array_FAI_counter); + + pqueue_bsf result = *pq_bsf; + free(pq_bsf); + return result; +} + +pqueue_bsf Fresh::FRESH_search_topk_DTW(ts_type *ts, node_list *nodelist, idx_t k) +{ + isax_index *index = this->index; + int warpWind = this->warping_window; + float *rawfile = this->database; + + ts_type *paa = (ts_type *)malloc(sizeof(ts_type) * index->settings->paa_segments); + paa_from_ts(ts, paa, index->settings->paa_segments, index->settings->ts_values_per_paa_segment); + + pqueue_bsf *pq_bsf = pqueue_bsf_init(k); + approximate_DTWtopk_inmemory(ts, paa, index, warpWind, pq_bsf, rawfile); + this->minimum_distance = pq_bsf->knn[k - 1]; + + ts_type *upperLemire = (ts_type *)malloc(sizeof(ts_type) * index->settings->timeseries_size); + ts_type *lowerLemire = (ts_type *)malloc(sizeof(ts_type) * index->settings->timeseries_size); + ts_type *paaU = (ts_type *)malloc(sizeof(ts_type) * index->settings->paa_segments); + ts_type *paaL = (ts_type *)malloc(sizeof(ts_type) * index->settings->paa_segments); + + lower_upper_lemire(ts, index->settings->timeseries_size, warpWind, lowerLemire, upperLemire); + paa_from_ts(upperLemire, paaU, index->settings->paa_segments, index->settings->ts_values_per_paa_segment); + paa_from_ts(lowerLemire, paaL, index->settings->paa_segments, index->settings->ts_values_per_paa_segment); + + int n_queues = (this->search_workers == 1) ? 1 : this->search_workers / 2; + int query_id = __sync_add_and_fetch(&g_fresh_query_id, 1); + + fresh_array_list_t *array_lists = (fresh_array_list_t *)malloc(sizeof(fresh_array_list_t) * n_queues); + for (int i = 0; i < n_queues; i++) + init_array_list_fresh(&array_lists[i], index->settings->root_nodes_size / n_queues); + + fresh_sorted_array_t *volatile *sorted_arrays = (fresh_sorted_array_t *volatile *)calloc(n_queues, sizeof(fresh_sorted_array_t *)); + volatile unsigned long *next_queue_data_pos = (volatile unsigned long *)calloc(n_queues, sizeof(unsigned long)); + volatile unsigned char *queue_finished = (volatile unsigned char *)calloc(n_queues, sizeof(unsigned char)); + volatile unsigned char *helper_queue_exist = (volatile unsigned char *)calloc(n_queues, sizeof(unsigned char)); + volatile unsigned long *sorted_array_FAI_counter = (volatile unsigned long *)calloc(n_queues, sizeof(unsigned long)); + + volatile int node_counter = 0; + + pthread_rwlock_t lock_bsf = PTHREAD_RWLOCK_INITIALIZER; + pthread_t threadid[this->search_workers]; + FRESH_workerdata workerdata[this->search_workers]; + + for (int i = 0; i < this->search_workers; i++) + { + workerdata[i].paa = paa; + workerdata[i].paaU = paaU; + workerdata[i].paaL = paaL; + workerdata[i].ts = ts; + workerdata[i].uo = upperLemire; + workerdata[i].lo = lowerLemire; + workerdata[i].index = index; + workerdata[i].minimum_distance = this->minimum_distance; + workerdata[i].pq_bsf = pq_bsf; + workerdata[i].lock_bsf = &lock_bsf; + workerdata[i].nodelist = nodelist->nlist; + workerdata[i].amountnode = nodelist->node_amount; + workerdata[i].node_counter = &node_counter; + workerdata[i].array_lists = array_lists; + workerdata[i].sorted_arrays = sorted_arrays; + workerdata[i].next_queue_data_pos = next_queue_data_pos; + workerdata[i].queue_finished = queue_finished; + workerdata[i].helper_queue_exist = helper_queue_exist; + workerdata[i].sorted_array_FAI_counter = sorted_array_FAI_counter; + workerdata[i].sorted_array_counter = nullptr; + workerdata[i].fai_queue_counters = nullptr; + workerdata[i].queue_bsf_distance = nullptr; + workerdata[i].workernumber = i; + workerdata[i].n_queues = n_queues; + workerdata[i].warpWind = warpWind; + workerdata[i].query_id = query_id; + workerdata[i].rawfile = rawfile; + } + + for (int i = 0; i < this->search_workers; i++) + pthread_create(&threadid[i], NULL, FRESH_topk_search_worker_DTW, (void *)&workerdata[i]); + for (int i = 0; i < this->search_workers; i++) + pthread_join(threadid[i], NULL); + + this->minimum_distance = pq_bsf->knn[k - 1]; + + for (int i = 0; i < n_queues; i++) + { + if (sorted_arrays[i]) + { + free(sorted_arrays[i]->data); + free(sorted_arrays[i]); + } + fresh_array_list_node_t *node = array_lists[i].Top; + while (node) + { + fresh_array_list_node_t *next = node->next; + free(node->data); + free(node); + node = next; + } + } + free(array_lists); + free((void *)sorted_arrays); + free((void *)next_queue_data_pos); + free((void *)queue_finished); + free((void *)helper_queue_exist); + free((void *)sorted_array_FAI_counter); + free(upperLemire); + free(lowerLemire); + free(paaU); + free(paaL); + free(paa); + + pqueue_bsf result = *pq_bsf; + free(pq_bsf); + return result; +} + +void Fresh::searchIndexL2Squared(const float *query, idx_t n_query, idx_t k, idx_t *I, float *D) +{ + ts_type *paa = (ts_type *)malloc(sizeof(ts_type) * index->settings->paa_segments); + + node_list nodelist; + nodelist.nlist = (isax_node **)malloc(sizeof(isax_node *) * pow(2, index->settings->paa_segments)); + nodelist.node_amount = 0; + isax_node *current_root_node = index->first_node; + while (current_root_node != NULL) + { + nodelist.nlist[nodelist.node_amount] = current_root_node; + current_root_node = current_root_node->next; + nodelist.node_amount++; + } + + for (idx_t q_loaded = 0; q_loaded < n_query; q_loaded++) + { + const float *ts = query + q_loaded * this->dim; + paa_from_ts(ts, paa, index->settings->paa_segments, index->settings->ts_values_per_paa_segment); + + pqueue_bsf result = FRESH_search_topk_L2Squared((float *)ts, paa, &nodelist, k); + + std::vector> pairs; + pairs.reserve(static_cast(k)); + for (idx_t ik = 0; ik < k; ik++) + { + if (result.position[ik] >= 0 && result.knn[ik] < FLT_MAX * 0.99f) + pairs.emplace_back(result.knn[ik], result.position[ik]); + } + std::sort(pairs.begin(), pairs.end(), [](const auto &a, const auto &b) { + if (a.first != b.first) return a.first < b.first; + return a.second < b.second; + }); + std::unordered_set seen_pos; + std::vector> uniq; + uniq.reserve(pairs.size()); + for (const auto &p : pairs) + { + if (seen_pos.insert(p.second).second) + uniq.push_back(p); + } + long last_pos = 0; + float last_dist = 0.0f; + for (idx_t ik = 0; ik < k; ik++) + { + if (ik < static_cast(uniq.size())) + { + last_dist = uniq[static_cast(ik)].first; + last_pos = uniq[static_cast(ik)].second; + } + I[q_loaded * k + ik] = static_cast(last_pos >= 0 ? last_pos : 0); + D[q_loaded * k + ik] = last_dist; + } + + free(result.position); + free(result.knn); + free(result.node); + } + + free(nodelist.nlist); + free(paa); + fprintf(stderr, ">>> Finished querying.\n"); + fflush(stdout); +} + +void Fresh::searchIndexDTW(const float *query, idx_t n_query, idx_t k, idx_t *I, float *D) +{ + isax_index *index = this->index; + + node_list nodelist; + nodelist.nlist = (isax_node **)malloc(sizeof(isax_node *) * pow(2, index->settings->paa_segments)); + nodelist.node_amount = 0; + isax_node *current_root_node = index->first_node; + while (current_root_node != NULL) + { + nodelist.nlist[nodelist.node_amount] = current_root_node; + current_root_node = current_root_node->next; + nodelist.node_amount++; + } + + for (idx_t q_loaded = 0; q_loaded < n_query; q_loaded++) + { + float *ts = (float *)&query[q_loaded * this->dim]; + pqueue_bsf result = FRESH_search_topk_DTW((float *)ts, &nodelist, k); + + std::vector> pairs; + pairs.reserve(static_cast(k)); + for (idx_t ik = 0; ik < k; ik++) + { + if (result.position[ik] >= 0 && result.knn[ik] < FLT_MAX * 0.99f) + pairs.emplace_back(result.knn[ik], result.position[ik]); + } + std::sort(pairs.begin(), pairs.end(), [](const auto &a, const auto &b) { + if (a.first != b.first) return a.first < b.first; + return a.second < b.second; + }); + std::unordered_set seen_pos; + std::vector> uniq; + uniq.reserve(pairs.size()); + for (const auto &p : pairs) + { + if (seen_pos.insert(p.second).second) + uniq.push_back(p); + } + long last_pos = 0; + float last_dist = 0.0f; + for (idx_t ik = 0; ik < k; ik++) + { + if (ik < static_cast(uniq.size())) + { + last_dist = uniq[static_cast(ik)].first; + last_pos = uniq[static_cast(ik)].second; + } + I[q_loaded * k + ik] = static_cast(last_pos >= 0 ? last_pos : 0); + D[q_loaded * k + ik] = last_dist; + } + + free(result.position); + free(result.knn); + free(result.node); + } + + free(nodelist.nlist); + fprintf(stderr, ">>> Finished querying.\n"); +} + +void Fresh::searchIndex(const float *query, idx_t n_query, idx_t k, idx_t *I, float *D) +{ + if (this->distance_type == DistanceType::L2_SQUARED) + searchIndexL2Squared(query, n_query, k, I, D); + else if (this->distance_type == DistanceType::DTW) + searchIndexDTW(query, n_query, k, I, D); + else + { + fprintf(stderr, "Error: Unsupported distance type for Fresh index.\n"); + exit(1); + } +} + +Fresh::~Fresh() +{ + if (owns_database && database != nullptr) + delete[] database; + + if (index != nullptr) + { + if (index->sax_cache != nullptr) + free(index->sax_cache); + if (index->answer != nullptr) + free(index->answer); + if (index->fbl != nullptr) + { + parallel_first_buffer_layer_lf *fbl = (parallel_first_buffer_layer_lf *)index->fbl; + for (int i = 0; i < fbl->number_of_buffers; i++) + { + parallel_fbl_soft_buffer_lf *sb = &fbl->soft_buffers[i]; + if (!sb->initialized) continue; + for (int k = 0; k < fbl->initial_buffer_size; k++) + { + if (sb->sax_records && sb->sax_records[k]) free(sb->sax_records[k]); + if (sb->pos_records && sb->pos_records[k]) free(sb->pos_records[k]); + } + free(sb->max_buffer_size); + free(sb->buffer_size); + free(sb->sax_records); + free(sb->pos_records); + } + free(fbl->soft_buffers); + free(fbl); + } + if (index->sax_file != nullptr) + fclose(index->sax_file); + free(index); + } + + if (index_settings != nullptr) + { + if (index_settings->bit_masks != nullptr) + free(index_settings->bit_masks); + if (index_settings->max_sax_cardinalities != nullptr) + free(index_settings->max_sax_cardinalities); + free(index_settings); + } +} + +} // namespace daisy + + diff --git a/lib/algos/Fresh.hpp b/lib/algos/Fresh.hpp new file mode 100644 index 0000000..fe2e1bf --- /dev/null +++ b/lib/algos/Fresh.hpp @@ -0,0 +1,175 @@ +#ifndef FRESH_HPP +#define FRESH_HPP + +#include "SimilaritySearchAlgorithm.hpp" +#include +#include +#include "../isax/iSAXIndex.hpp" +#include "../isax/iSAXPqueue.hpp" + +#define FRESH_CASPTR(A,B,C) __sync_bool_compare_and_swap((long *)(A), (long)(B), (long)(C)) +#define FRESH_CASULONG(A,B,C) __sync_bool_compare_and_swap((unsigned long *)(A), (unsigned long)(B), (unsigned long)(C)) + +namespace daisy +{ + + typedef struct fresh_array_element + { + float distance; + isax_node *node; + } fresh_array_element_t; + + typedef struct fresh_array_list_node + { + fresh_array_element_t *data; + int num_node; + struct fresh_array_list_node *next; + } fresh_array_list_node_t; + + typedef struct fresh_array_list + { + fresh_array_list_node_t *Top; + int element_size; + } fresh_array_list_t; + + typedef struct fresh_sorted_array + { + fresh_array_element_t *data; + int num_elements; + } fresh_sorted_array_t; + + typedef struct parallel_fbl_soft_buffer_lf + { + isax_node *volatile node; + sax_type **sax_records; + file_position_type **pos_records; + volatile unsigned char initialized; + int *max_buffer_size; + int *buffer_size; + volatile unsigned char processed; + volatile unsigned long next_iSAX_group; + root_mask_type mask; + volatile unsigned char **iSAX_processed; + volatile unsigned char recBuf_helpers_exist; + } parallel_fbl_soft_buffer_lf; + + typedef struct parallel_first_buffer_layer_lf + { + int number_of_buffers; + int initial_buffer_size; + int max_total_size; + int current_record_index; + parallel_fbl_soft_buffer_lf *soft_buffers; + } parallel_first_buffer_layer_lf; + + struct FreshConfig + { + int search_workers = 4; + int index_workers = 2; + int warping_window = 10; + int leaf_size = 2000; + int paa_segments = 16; + }; + + // FreSH's lock-free workerdata. bsf_result_p replaced by pq_bsf + lock_bsf + // because DaiSy requires top-k; FreSH's CAS BSF update is 1-NN only. + typedef struct FRESH_workerdata + { + ts_type *paa, *paaU, *paaL, *ts, *uo, *lo; + fresh_array_list_t *array_lists; + fresh_sorted_array_t *volatile *sorted_arrays; + volatile unsigned char *queue_finished; + volatile unsigned char *helper_queue_exist; + volatile int *fai_queue_counters; + volatile float **queue_bsf_distance; + isax_index *index; + float minimum_distance; + pqueue_bsf *pq_bsf; + pthread_rwlock_t *lock_bsf; + volatile int *node_counter; + volatile unsigned long *sorted_array_counter; + volatile unsigned long *sorted_array_FAI_counter; + isax_node **nodelist; + int amountnode; + int workernumber; + int n_queues; + int warpWind; + volatile unsigned long *next_queue_data_pos; + int query_id; + pthread_barrier_t *wait_tree_pruning_phase_to_finish; + pthread_barrier_t *wait_process_queue; + float *rawfile; + } FRESH_workerdata; + + void *FRESH_topk_search_worker_L2Squared(void *rfdata); + void *FRESH_topk_search_worker_DTW(void *rfdata); + + class Fresh : public SimilaritySearchAlgorithm + { + private: + int read_block_length = 100000; + int search_workers = 4; + int index_workers = 2; + bool owns_database = false; + + pqueue_bsf FRESH_search_topk_L2Squared(ts_type *ts, ts_type *paa, node_list *nodelist, idx_t k); + pqueue_bsf FRESH_search_topk_DTW(ts_type *ts, node_list *nodelist, idx_t k); + + void searchIndexL2Squared(const float *query, idx_t n_query, idx_t k, idx_t *I, float *D); + void searchIndexDTW(const float *query, idx_t n_query, idx_t k, idx_t *I, float *D); + + public: + Fresh(DistanceType distance_type); + Fresh(DistanceType distance_type, const FreshConfig &config); + + void setWarpingWindow(int w) { warping_window = w; } + void setWarpWindow(int w) { warping_window = w; } + + 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("Fresh requires in-memory data. Use buildIndex(database, n_database, dim) instead."); + } + + void searchIndex(const float *query, idx_t n_query, idx_t k, idx_t *I, float *D) override; + + int getNumThreads() const { return SimilaritySearchAlgorithm::num_threads; } + void setNumThreads(int n) { + SimilaritySearchAlgorithm::num_threads = n; + search_workers = n; + } + int getPaaSegments() const { return paa_segments; } + void setPaaSegments(int n) { paa_segments = n; } + int getSaxCardinality() const { return sax_cardinality; } + void setSaxCardinality(int n) { sax_cardinality = n; } + int getLeafSize() const { return leaf_size; } + void setLeafSize(int n) { leaf_size = n; } + int getMinLeafSize() const { return min_leaf_size; } + void setMinLeafSize(int n) { min_leaf_size = n; } + int getInitialLblSize() const { return initial_lbl_size; } + void setInitialLblSize(int n) { initial_lbl_size = n; } + int getFlushLimit() const { return flush_limit; } + void setFlushLimit(int n) { flush_limit = n; } + int getInitialFblSize() const { return initial_fbl_size; } + void setInitialFblSize(int n) { initial_fbl_size = n; } + int getTotalLoadedLeaves() const { return total_loaded_leaves; } + void setTotalLoadedLeaves(int n) { total_loaded_leaves = n; } + int getTightBound() const { return tight_bound; } + void setTightBound(int n) { tight_bound = n; } + int getSearchWorkers() const { return search_workers; } + void setSearchWorkers(int n) { search_workers = n; } + int getIndexWorkers() const { return index_workers; } + void setIndexWorkers(int n) { index_workers = n; } + int getReadBlockLength() const { return read_block_length; } + void setReadBlockLength(int n) { read_block_length = n; } + int getWarpingWindow() const { return warping_window; } + + ~Fresh(); + }; + +} + +#endif diff --git a/lib/daisy.hpp b/lib/daisy.hpp index d34ac2c..b4da494 100644 --- a/lib/daisy.hpp +++ b/lib/daisy.hpp @@ -10,5 +10,6 @@ #include "algos/Sofa.hpp" #include "algos/Hercules.hpp" #include "algos/DumpyOS.hpp" +#include "algos/Fresh.hpp" #endif diff --git a/lib/isax/iSAXIndex.cpp b/lib/isax/iSAXIndex.cpp index 0b996a5..9d97efd 100644 --- a/lib/isax/iSAXIndex.cpp +++ b/lib/isax/iSAXIndex.cpp @@ -471,6 +471,7 @@ namespace daisy node->buffer = init_node_buffer(initial_buffer_size); node->mask = 0; node->wedges = NULL; + node->processed = 0; return node; } diff --git a/lib/isax/iSAXIndex.hpp b/lib/isax/iSAXIndex.hpp index 87615be..ad6adb9 100644 --- a/lib/isax/iSAXIndex.hpp +++ b/lib/isax/iSAXIndex.hpp @@ -140,6 +140,9 @@ namespace daisy // Wedges ts_type *wedges; + // FreSH-specific: deduplicates subtree visits during lock-free pruning, unused by other algorithms + volatile unsigned long processed; + } isax_node; typedef struct fbl_soft_buffer diff --git a/pybinds/setup.cpp b/pybinds/setup.cpp index cb56223..be35e1c 100644 --- a/pybinds/setup.cpp +++ b/pybinds/setup.cpp @@ -25,6 +25,7 @@ #include "../lib/algos/DataSource.hpp" #include "../lib/algos/Hercules.hpp" #include "../lib/algos/DumpyOS.hpp" +#include "../lib/algos/Fresh.hpp" #ifdef SOFA_FFTW_ENABLED #if SOFA_FFTW_ENABLED != 0 #include "../lib/algos/Sofa.hpp" @@ -638,4 +639,74 @@ PYBIND11_MODULE(_core, m) pybind11::array_t({n_query, k}, indices.data()), pybind11::array_t({n_query, k}, distances.data()) ); }, "Search the DumpyOS index and return (indices, distances)"); + + ////// FRESH ////// + pybind11::class_(m, "Fresh", "FreSH lock-free iSAX-based time series similarity index (SRDS 2023)") + .def(pybind11::init(), "Create a new Fresh instance with the given distance metric") + + // Getters + .def("getNumThreads", &daisy::Fresh::getNumThreads, "Get the number of search threads") + .def("getPaaSegments", &daisy::Fresh::getPaaSegments, "Get the number of PAA segments used in SAX transformation") + .def("getSaxCardinality", &daisy::Fresh::getSaxCardinality, "Get the cardinality of SAX symbols") + .def("getLeafSize", &daisy::Fresh::getLeafSize, "Get the maximum leaf size in the index tree") + .def("getMinLeafSize", &daisy::Fresh::getMinLeafSize, "Get the minimum number of entries per leaf") + .def("getInitialLblSize", &daisy::Fresh::getInitialLblSize, "Get the initial size of the lower-bound buffer") + .def("getFlushLimit", &daisy::Fresh::getFlushLimit, "Get the flush limit before writing to disk") + .def("getInitialFblSize", &daisy::Fresh::getInitialFblSize, "Get the initial full-buffer size") + .def("getTotalLoadedLeaves", &daisy::Fresh::getTotalLoadedLeaves, "Get the total number of leaves loaded") + .def("getTightBound", &daisy::Fresh::getTightBound, "Check whether tight bounds are enabled") + .def("getSearchWorkers", &daisy::Fresh::getSearchWorkers, "Get number of worker threads used for search") + .def("getIndexWorkers", &daisy::Fresh::getIndexWorkers, "Get number of worker threads used for indexing") + .def("getReadBlockLength", &daisy::Fresh::getReadBlockLength, "Get block size for reading the time series data") + .def("getWarpingWindow", &daisy::Fresh::getWarpingWindow, "Get the DTW warping window constraint") + + // Setters + .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("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") + .def("setFlushLimit", &daisy::Fresh::setFlushLimit, "Set the flush limit") + .def("setInitialFblSize", &daisy::Fresh::setInitialFblSize, "Set the initial FBL size") + .def("setTotalLoadedLeaves", &daisy::Fresh::setTotalLoadedLeaves, "Set the number of total loaded leaves") + .def("setTightBound", &daisy::Fresh::setTightBound, "Enable or disable tight bounds") + .def("setSearchWorkers", &daisy::Fresh::setSearchWorkers, "Set the number of worker threads for search") + .def("setIndexWorkers", &daisy::Fresh::setIndexWorkers, "Set the number of worker threads for indexing") + .def("setReadBlockLength", &daisy::Fresh::setReadBlockLength, "Set the length of each read block") + .def("setWarpingWindow", &daisy::Fresh::setWarpingWindow, "Set the warping window size for DTW") + + // Build the index from a 2D NumPy array + .def("buildIndex", [](daisy::Fresh &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 Fresh index from a 2D float32 NumPy array") + + // Search the index with query array and return top-k results + .def("searchIndex", [](daisy::Fresh &self, pybind11::array_t query, daisy::idx_t k) + { + pybind11::buffer_info query_buf = query.request(); + if (query_buf.ndim != 2) + throw std::runtime_error("Query array must be 2D"); + if (k <= 0) + throw std::runtime_error("k must be positive"); + + const daisy::idx_t n_query = query_buf.shape[0]; + + std::vector indices(n_query * k); + std::vector distances(n_query * k); + + self.searchIndex(static_cast(query_buf.ptr), n_query, k, indices.data(), distances.data()); + + return pybind11::make_tuple( + pybind11::array_t({n_query, k}, indices.data()), + pybind11::array_t({n_query, k}, distances.data()) + ); }, "Search the Fresh index using queries and return (indices, distances)"); } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8ad95bd..0c353c6 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -678,4 +678,31 @@ target_include_directories(test_DumpyOS_DTW gtest_discover_tests( test_DumpyOS_DTW WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} +) + +# ////// FRESH L2Square ////// +add_executable( + test_Fresh_L2Square + test_Fresh_L2Square.cpp + test_utils.cpp +) + +target_link_libraries( + test_Fresh_L2Square + PRIVATE + GTest::gtest_main + dino_lib + commons_lib + stdc++fs +) + +target_include_directories(test_Fresh_L2Square + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons +) + +gtest_discover_tests( + test_Fresh_L2Square + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} ) \ No newline at end of file diff --git a/tests/test_Fresh_L2Square.cpp b/tests/test_Fresh_L2Square.cpp new file mode 100644 index 0000000..aabb0f8 --- /dev/null +++ b/tests/test_Fresh_L2Square.cpp @@ -0,0 +1,45 @@ +#include "test_utils.hpp" +#include "../commons/test_bm_utils.hpp" +#include "../commons/paramSetup.hpp" + +std::string prefix = "bruteForce"; + +TEST_P(FreshParameterizedTest, AllConfigurations) +{ + const SSTestConfig &config = GetParam(); + daisy::DistanceType dist_L2Squared = daisy::DistanceType::L2_SQUARED; + for (int i = 0; i < 3; ++i) + { + daisy::Fresh search(dist_L2Squared); + + std::string gt_I_path = config.gt_I_prefix + std::to_string(config.k_value) + ".txt"; + std::string gt_D_path = config.gt_D_prefix + std::to_string(config.k_value) + ".txt"; + + runSST( + &search, + prefix, + gt_I_path, + gt_D_path, + config.dataset_path, + config.query_path, + config.thread_count); + } + +} + +INSTANTIATE_TEST_SUITE_P( + FreshTests, + FreshParameterizedTest, + ::testing::ValuesIn(test_configs), + [](const ::testing::TestParamInfo &info) + { + return info.param.name + "_k" + std::to_string(info.param.k_value) + + "_thread" + std::to_string(info.param.thread_count) + + "_idx" + std::to_string(info.index); + }); + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/tests/test_utils.hpp b/tests/test_utils.hpp index 8ac8a56..9276e1c 100644 --- a/tests/test_utils.hpp +++ b/tests/test_utils.hpp @@ -16,6 +16,7 @@ #include "../lib/algos/Sofa.hpp" #include "../lib/algos/Hercules.hpp" #include "../lib/algos/DumpyOS.hpp" +#include "../lib/algos/Fresh.hpp" class SimilaritySearchTest : public ::testing::Test { @@ -163,4 +164,14 @@ class DumpyOSDTWParameterizedTest : public SimilaritySearchTest, static void TearDownTestSuite() {} }; +class FreshParameterizedTest : public SimilaritySearchTest, + public ::testing::WithParamInterface +{ +protected: + using SimilaritySearchTest::runSST; + + static void SetUpTestSuite() {} + static void TearDownTestSuite() {} +}; + #endif From 8d8ceb851d38dda74a97d652c24af7ec270ef638 Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Wed, 1 Jul 2026 15:32:02 +0200 Subject: [PATCH 11/19] fix bug in fresh --- lib/algos/Fresh.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/algos/Fresh.cpp b/lib/algos/Fresh.cpp index 3e10419..44f4d62 100644 --- a/lib/algos/Fresh.cpp +++ b/lib/algos/Fresh.cpp @@ -1274,10 +1274,10 @@ Fresh::~Fresh() { parallel_fbl_soft_buffer_lf *sb = &fbl->soft_buffers[i]; if (!sb->initialized) continue; - for (int k = 0; k < fbl->initial_buffer_size; k++) + for (int w = 0; w < this->index_workers; w++) { - if (sb->sax_records && sb->sax_records[k]) free(sb->sax_records[k]); - if (sb->pos_records && sb->pos_records[k]) free(sb->pos_records[k]); + if (sb->sax_records && sb->sax_records[w]) free(sb->sax_records[w]); + if (sb->pos_records && sb->pos_records[w]) free(sb->pos_records[w]); } free(sb->max_buffer_size); free(sb->buffer_size); From 1e6ad6fdd8c85c335ebfa79089ead5621ce86b70 Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Wed, 1 Jul 2026 21:32:32 +0200 Subject: [PATCH 12/19] fix fresh: set is_leaf=1 on node init, revert isax changes --- lib/algos/Fresh.cpp | 22 +++++++--------------- lib/isax/iSAXIndex.cpp | 1 - lib/isax/iSAXIndex.hpp | 3 --- 3 files changed, 7 insertions(+), 19 deletions(-) diff --git a/lib/algos/Fresh.cpp b/lib/algos/Fresh.cpp index 44f4d62..87806b8 100644 --- a/lib/algos/Fresh.cpp +++ b/lib/algos/Fresh.cpp @@ -391,6 +391,7 @@ static void *indexCreationWorkerFresh(void *transferdata) if (!current_fbl_node->node) { isax_node *root = isax_root_node_init(current_fbl_node->mask, index->settings->initial_leaf_buffer_size); + root->is_leaf = 1; if (!FRESH_CASPTR(¤t_fbl_node->node, NULL, root)) free(root); } @@ -461,7 +462,7 @@ static void add_to_array_data_lf(float *paa, isax_node *node, isax_index *index, volatile unsigned long *next_queue_data_pos, int n_queues, const int query_id) { - if (node->processed == (unsigned long)query_id) + if (node->isax_values == NULL) return; float distance = minidist_paa_to_isax(paa, node->isax_values, @@ -472,7 +473,7 @@ static void add_to_array_data_lf(float *paa, isax_node *node, isax_index *index, MINVAL, MAXVAL, index->settings->mindist_sqrt); - if (distance < bsf && node->processed < (unsigned long)query_id) + if (distance < bsf) { if (node->is_leaf) { @@ -484,15 +485,12 @@ static void add_to_array_data_lf(float *paa, isax_node *node, isax_index *index, } else { - if (node->left_child->isax_cardinalities != NULL && node->left_child->processed < (unsigned long)query_id) + if (node->left_child->isax_cardinalities != NULL) add_to_array_data_lf(paa, node->left_child, index, bsf, array_lists, tnumber, next_queue_data_pos, n_queues, query_id); - if (node->right_child->isax_cardinalities != NULL && node->right_child->processed < (unsigned long)query_id) + if (node->right_child->isax_cardinalities != NULL) add_to_array_data_lf(paa, node->right_child, index, bsf, array_lists, tnumber, next_queue_data_pos, n_queues, query_id); } } - - if (node->processed < (unsigned long)query_id) - node->processed = query_id; } static int compare_sorted_array_items(const void *a, const void *b) @@ -642,10 +640,7 @@ void *FRESH_topk_search_worker_L2Squared(void *rfdata) // A.1. help unprocessed subtrees for (int i = 0; i < wd->amountnode; i++) { - isax_node *node = wd->nodelist[i]; - if (node->processed == (unsigned long)query_id) - continue; - add_to_array_data_lf(paa, node, index, bsfdistance, + add_to_array_data_lf(paa, wd->nodelist[i], index, bsfdistance, wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); } @@ -701,10 +696,7 @@ void *FRESH_topk_search_worker_DTW(void *rfdata) // A.1. help unprocessed subtrees for (int i = 0; i < wd->amountnode; i++) { - isax_node *node = wd->nodelist[i]; - if (node->processed == (unsigned long)query_id) - continue; - add_to_array_data_lf(paa, node, index, bsfdistance, + add_to_array_data_lf(paa, wd->nodelist[i], index, bsfdistance, wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); } diff --git a/lib/isax/iSAXIndex.cpp b/lib/isax/iSAXIndex.cpp index 9d97efd..0b996a5 100644 --- a/lib/isax/iSAXIndex.cpp +++ b/lib/isax/iSAXIndex.cpp @@ -471,7 +471,6 @@ namespace daisy node->buffer = init_node_buffer(initial_buffer_size); node->mask = 0; node->wedges = NULL; - node->processed = 0; return node; } diff --git a/lib/isax/iSAXIndex.hpp b/lib/isax/iSAXIndex.hpp index ad6adb9..87615be 100644 --- a/lib/isax/iSAXIndex.hpp +++ b/lib/isax/iSAXIndex.hpp @@ -140,9 +140,6 @@ namespace daisy // Wedges ts_type *wedges; - // FreSH-specific: deduplicates subtree visits during lock-free pruning, unused by other algorithms - volatile unsigned long processed; - } isax_node; typedef struct fbl_soft_buffer From 70215819572c2af65a7098a992e25ee0fd55597c Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Wed, 1 Jul 2026 21:54:32 +0200 Subject: [PATCH 13/19] fix fresh search: fix lf struct layout, build first_node list after indexing --- lib/algos/Fresh.cpp | 19 +++++++++++++++++++ lib/algos/Fresh.hpp | 14 +++++++++----- 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/lib/algos/Fresh.cpp b/lib/algos/Fresh.cpp index 87806b8..ff89dc9 100644 --- a/lib/algos/Fresh.cpp +++ b/lib/algos/Fresh.cpp @@ -68,6 +68,8 @@ static parallel_first_buffer_layer_lf *initialize_pRecBuf_lf(int initial_buffer_ fbl->initial_buffer_size = initial_buffer_size; fbl->max_total_size = max_total_size; fbl->current_record_index = 0; + fbl->current_record = NULL; + fbl->hard_buffer = NULL; fbl->soft_buffers = (parallel_fbl_soft_buffer_lf *)calloc(number_of_buffers, sizeof(parallel_fbl_soft_buffer_lf)); return fbl; } @@ -880,6 +882,23 @@ void Fresh::buildIndex(DataSource *data_source) for (int i = 0; i < this->index_workers; i++) pthread_join(threadid[i], NULL); + // Build index->first_node linked list from FBL so that search functions can traverse all roots. + { + parallel_first_buffer_layer_lf *fbl = (parallel_first_buffer_layer_lf *)index->fbl; + for (int i = 0; i < fbl->number_of_buffers; i++) + { + if (!fbl->soft_buffers[i].initialized || !fbl->soft_buffers[i].node) + continue; + isax_node *node = fbl->soft_buffers[i].node; + node->next = index->first_node; + node->previous = NULL; + if (index->first_node) + index->first_node->previous = node; + index->first_node = node; + __sync_fetch_and_add(&(index->root_nodes), 1); + } + } + __sync_fetch_and_add(&(index->total_records), this->n_database); index->sax_cache_size = index->total_records; fprintf(stderr, ">>> Finished indexing\n"); diff --git a/lib/algos/Fresh.hpp b/lib/algos/Fresh.hpp index fe2e1bf..f1e264f 100644 --- a/lib/algos/Fresh.hpp +++ b/lib/algos/Fresh.hpp @@ -38,27 +38,31 @@ namespace daisy int num_elements; } fresh_sorted_array_t; + // Must stay layout-compatible with parallel_fbl_soft_buffer (same field offsets, same 56-byte size) + // because iSAXSearch.cpp casts index->fbl through parallel_first_buffer_layer * and accesses + // .initialized (offset 24) and .node (offset 0) via the MESSI slot type. typedef struct parallel_fbl_soft_buffer_lf { isax_node *volatile node; sax_type **sax_records; file_position_type **pos_records; - volatile unsigned char initialized; + int initialized; + int _pad; int *max_buffer_size; int *buffer_size; - volatile unsigned char processed; - volatile unsigned long next_iSAX_group; root_mask_type mask; - volatile unsigned char **iSAX_processed; - volatile unsigned char recBuf_helpers_exist; } parallel_fbl_soft_buffer_lf; + // Must stay layout-compatible with parallel_first_buffer_layer so that + // iSAXSearch.cpp's (parallel_first_buffer_layer *)(index->fbl)->soft_buffers lands at offset 32. typedef struct parallel_first_buffer_layer_lf { int number_of_buffers; int initial_buffer_size; int max_total_size; int current_record_index; + char *current_record; + char *hard_buffer; parallel_fbl_soft_buffer_lf *soft_buffers; } parallel_first_buffer_layer_lf; From 13508e44ac2a3bbfac921f781d959ddabb0adeb6 Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Wed, 1 Jul 2026 22:36:46 +0200 Subject: [PATCH 14/19] fix fresh bug --- lib/algos/Fresh.cpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/lib/algos/Fresh.cpp b/lib/algos/Fresh.cpp index ff89dc9..92ff025 100644 --- a/lib/algos/Fresh.cpp +++ b/lib/algos/Fresh.cpp @@ -646,6 +646,11 @@ void *FRESH_topk_search_worker_L2Squared(void *rfdata) wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); } + // All workers must finish populating queues before any worker sorts them. + // Without this barrier a fast worker can snapshot next_queue_data_pos while + // slower workers are still appending, silently dropping those entries. + pthread_barrier_wait(wd->wait_tree_pruning_phase_to_finish); + // A.2. create my sorted array int my_pq = wd->workernumber % wd->n_queues; if (!wd->sorted_arrays[my_pq] && wd->next_queue_data_pos[my_pq]) @@ -702,6 +707,8 @@ void *FRESH_topk_search_worker_DTW(void *rfdata) wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); } + pthread_barrier_wait(wd->wait_tree_pruning_phase_to_finish); + // A.2. create my sorted array int my_pq = wd->workernumber % wd->n_queues; if (!wd->sorted_arrays[my_pq] && wd->next_queue_data_pos[my_pq]) @@ -951,6 +958,9 @@ pqueue_bsf Fresh::FRESH_search_topk_L2Squared(ts_type *ts, ts_type *paa, node_li volatile int node_counter = 0; + pthread_barrier_t wait_tree_pruning_phase_to_finish; + pthread_barrier_init(&wait_tree_pruning_phase_to_finish, NULL, this->search_workers); + pthread_rwlock_t lock_bsf = PTHREAD_RWLOCK_INITIALIZER; pthread_t threadid[this->search_workers]; FRESH_workerdata workerdata[this->search_workers]; @@ -975,6 +985,7 @@ pqueue_bsf Fresh::FRESH_search_topk_L2Squared(ts_type *ts, ts_type *paa, node_li workerdata[i].sorted_array_counter = nullptr; workerdata[i].fai_queue_counters = nullptr; workerdata[i].queue_bsf_distance = nullptr; + workerdata[i].wait_tree_pruning_phase_to_finish = &wait_tree_pruning_phase_to_finish; workerdata[i].workernumber = i; workerdata[i].n_queues = n_queues; workerdata[i].query_id = query_id; @@ -986,6 +997,8 @@ pqueue_bsf Fresh::FRESH_search_topk_L2Squared(ts_type *ts, ts_type *paa, node_li for (int i = 0; i < this->search_workers; i++) pthread_join(threadid[i], NULL); + pthread_barrier_destroy(&wait_tree_pruning_phase_to_finish); + this->minimum_distance = pq_bsf->knn[k - 1]; for (int i = 0; i < n_queues; i++) @@ -1053,6 +1066,9 @@ pqueue_bsf Fresh::FRESH_search_topk_DTW(ts_type *ts, node_list *nodelist, idx_t volatile int node_counter = 0; + pthread_barrier_t wait_tree_pruning_phase_to_finish; + pthread_barrier_init(&wait_tree_pruning_phase_to_finish, NULL, this->search_workers); + pthread_rwlock_t lock_bsf = PTHREAD_RWLOCK_INITIALIZER; pthread_t threadid[this->search_workers]; FRESH_workerdata workerdata[this->search_workers]; @@ -1081,6 +1097,7 @@ pqueue_bsf Fresh::FRESH_search_topk_DTW(ts_type *ts, node_list *nodelist, idx_t workerdata[i].sorted_array_counter = nullptr; workerdata[i].fai_queue_counters = nullptr; workerdata[i].queue_bsf_distance = nullptr; + workerdata[i].wait_tree_pruning_phase_to_finish = &wait_tree_pruning_phase_to_finish; workerdata[i].workernumber = i; workerdata[i].n_queues = n_queues; workerdata[i].warpWind = warpWind; @@ -1093,6 +1110,8 @@ pqueue_bsf Fresh::FRESH_search_topk_DTW(ts_type *ts, node_list *nodelist, idx_t for (int i = 0; i < this->search_workers; i++) pthread_join(threadid[i], NULL); + pthread_barrier_destroy(&wait_tree_pruning_phase_to_finish); + this->minimum_distance = pq_bsf->knn[k - 1]; for (int i = 0; i < n_queues; i++) From c88ea50724c9ab60400e60c6fcb0c9e8be14191f Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Wed, 1 Jul 2026 22:49:10 +0200 Subject: [PATCH 15/19] Fix persisting bug: include minidist=0 leaf nodes in FRESH search sorted arrays --- lib/algos/Fresh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/algos/Fresh.cpp b/lib/algos/Fresh.cpp index 92ff025..ed30e04 100644 --- a/lib/algos/Fresh.cpp +++ b/lib/algos/Fresh.cpp @@ -522,7 +522,7 @@ static void create_sorted_array_from_data_queue(int pq_id, fresh_array_list_t *a array_elements_traversed < array_elements && !sorted_arrays[pq_id]; k++, array_elements_traversed++) { - if (bucket->data[k].node && bucket->data[k].distance > 0) + if (bucket->data[k].node && bucket->data[k].distance >= 0) { local_sa->data[j].node = bucket->data[k].node; local_sa->data[j].distance = bucket->data[k].distance; From 97af2a6e1f622751b63488aabe22480b9b752f84 Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Wed, 1 Jul 2026 23:05:20 +0200 Subject: [PATCH 16/19] Fix fresh sorted-array traversal order and remove duplication --- lib/algos/Fresh.cpp | 39 +++++++++++++++++---------------------- 1 file changed, 17 insertions(+), 22 deletions(-) diff --git a/lib/algos/Fresh.cpp b/lib/algos/Fresh.cpp index ed30e04..8f4f7d7 100644 --- a/lib/algos/Fresh.cpp +++ b/lib/algos/Fresh.cpp @@ -515,21 +515,33 @@ static void create_sorted_array_from_data_queue(int pq_id, fresh_array_list_t *a size_t j = 0; int array_buckets_num = array_lists[pq_id].Top->num_node + 1; unsigned long array_elements_traversed = 0; - fresh_array_list_node_t *bucket = array_lists[pq_id].Top; - for (int i = 0; i < array_buckets_num; i++, bucket = bucket->next) + + // The list is a stack: Top is the newest bucket (highest num_node), .next goes + // toward older buckets. Positions increase with num_node (bucket i holds positions + // i*element_size .. (i+1)*element_size-1). We must traverse oldest-first so the + // array_elements_traversed limit cuts off at the END of the newest bucket, not the + // middle of the oldest one. + fresh_array_list_node_t **buckets = (fresh_array_list_node_t **)malloc( + (size_t)array_buckets_num * sizeof(fresh_array_list_node_t *)); + fresh_array_list_node_t *b = array_lists[pq_id].Top; + for (int i = array_buckets_num - 1; i >= 0; i--, b = b->next) + buckets[i] = b; + + for (int i = 0; i < array_buckets_num; i++) { for (int k = 0; k < array_lists[pq_id].element_size && array_elements_traversed < array_elements && !sorted_arrays[pq_id]; k++, array_elements_traversed++) { - if (bucket->data[k].node && bucket->data[k].distance >= 0) + if (buckets[i]->data[k].node && buckets[i]->data[k].distance >= 0) { - local_sa->data[j].node = bucket->data[k].node; - local_sa->data[j].distance = bucket->data[k].distance; + local_sa->data[j].node = buckets[i]->data[k].node; + local_sa->data[j].distance = buckets[i]->data[k].distance; j++; } } } + free(buckets); local_sa->num_elements = j; if (!sorted_arrays[pq_id]) @@ -639,16 +651,6 @@ void *FRESH_topk_search_worker_L2Squared(void *rfdata) wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); } - // A.1. help unprocessed subtrees - for (int i = 0; i < wd->amountnode; i++) - { - add_to_array_data_lf(paa, wd->nodelist[i], index, bsfdistance, - wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); - } - - // All workers must finish populating queues before any worker sorts them. - // Without this barrier a fast worker can snapshot next_queue_data_pos while - // slower workers are still appending, silently dropping those entries. pthread_barrier_wait(wd->wait_tree_pruning_phase_to_finish); // A.2. create my sorted array @@ -700,13 +702,6 @@ void *FRESH_topk_search_worker_DTW(void *rfdata) wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); } - // A.1. help unprocessed subtrees - for (int i = 0; i < wd->amountnode; i++) - { - add_to_array_data_lf(paa, wd->nodelist[i], index, bsfdistance, - wd->array_lists, &tnumber, wd->next_queue_data_pos, wd->n_queues, query_id); - } - pthread_barrier_wait(wd->wait_tree_pruning_phase_to_finish); // A.2. create my sorted array From 2e5c269f835457017a3f1a889cf697f7034d2dfa Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Wed, 1 Jul 2026 23:10:35 +0200 Subject: [PATCH 17/19] Fix fresh sorted-array traversal order and remove duplication --- lib/algos/Fresh.cpp | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/lib/algos/Fresh.cpp b/lib/algos/Fresh.cpp index 8f4f7d7..e88a0b4 100644 --- a/lib/algos/Fresh.cpp +++ b/lib/algos/Fresh.cpp @@ -349,7 +349,6 @@ static void *indexCreationWorkerFresh(void *transferdata) sax_type *sax = (sax_type *)malloc(sizeof(sax_type) * index->settings->paa_segments); - // Phase 1: summarization — FAI on blocks, Refresh within each block unsigned long block_num; while (!*input_data->all_blocks_processed) { @@ -364,7 +363,6 @@ static void *indexCreationWorkerFresh(void *transferdata) pthread_barrier_wait(input_data->wait_summaries_to_compute); - // Phase 2: tree population — FAI on FBL slots, serial add_record_to_node_inmemory isax_node_record *r = (isax_node_record *)malloc(sizeof(isax_node_record)); while (1) @@ -411,7 +409,6 @@ static void *indexCreationWorkerFresh(void *transferdata) return NULL; } -// ── Phase 3: lock-free candidate queue (array list + sorted array) ──────────── static void init_array_list_fresh(fresh_array_list_t *list, int size) { @@ -554,7 +551,6 @@ static void create_sorted_array_from_data_queue(int pq_id, fresh_array_list_t *a } } -// ── Phase 4: refinement (help sorted arrays) ────────────────────────────────── static int process_sorted_array_element_L2(fresh_array_element_t *n, FRESH_workerdata *wd) { @@ -641,7 +637,6 @@ void *FRESH_topk_search_worker_L2Squared(void *rfdata) const int query_id = wd->query_id; int tnumber = rand() % wd->n_queues; - // A. populate arrays while (1) { int current_root_node_number = __sync_fetch_and_add(wd->node_counter, 1); @@ -653,24 +648,20 @@ void *FRESH_topk_search_worker_L2Squared(void *rfdata) pthread_barrier_wait(wd->wait_tree_pruning_phase_to_finish); - // A.2. create my sorted array int my_pq = wd->workernumber % wd->n_queues; if (!wd->sorted_arrays[my_pq] && wd->next_queue_data_pos[my_pq]) create_sorted_array_from_data_queue(my_pq, wd->array_lists, wd->next_queue_data_pos, wd->sorted_arrays); - // B. process my sorted array fresh_sorted_array_t *my_array = wd->sorted_arrays[my_pq]; if (my_array && !wd->queue_finished[my_pq]) help_sorted_array_L2(wd, my_array, my_pq); - // A.3. help create other sorted arrays for (int i = 0; i < wd->n_queues; i++) { if (!wd->sorted_arrays[i] && wd->next_queue_data_pos[i]) create_sorted_array_from_data_queue(i, wd->array_lists, wd->next_queue_data_pos, wd->sorted_arrays); } - // B.3. help process other sorted arrays for (int i = 0; i < wd->n_queues; i++) { if (wd->queue_finished[i] || !wd->sorted_arrays[i]) @@ -692,7 +683,6 @@ void *FRESH_topk_search_worker_DTW(void *rfdata) const int query_id = wd->query_id; int tnumber = rand() % wd->n_queues; - // A. populate arrays while (1) { int current_root_node_number = __sync_fetch_and_add(wd->node_counter, 1); @@ -704,24 +694,20 @@ void *FRESH_topk_search_worker_DTW(void *rfdata) pthread_barrier_wait(wd->wait_tree_pruning_phase_to_finish); - // A.2. create my sorted array int my_pq = wd->workernumber % wd->n_queues; if (!wd->sorted_arrays[my_pq] && wd->next_queue_data_pos[my_pq]) create_sorted_array_from_data_queue(my_pq, wd->array_lists, wd->next_queue_data_pos, wd->sorted_arrays); - // B. process my sorted array fresh_sorted_array_t *my_array = wd->sorted_arrays[my_pq]; if (my_array && !wd->queue_finished[my_pq]) help_sorted_array_DTW(wd, my_array, my_pq); - // A.3. help create other sorted arrays for (int i = 0; i < wd->n_queues; i++) { if (!wd->sorted_arrays[i] && wd->next_queue_data_pos[i]) create_sorted_array_from_data_queue(i, wd->array_lists, wd->next_queue_data_pos, wd->sorted_arrays); } - // B.3. help process other sorted arrays for (int i = 0; i < wd->n_queues; i++) { if (wd->queue_finished[i] || !wd->sorted_arrays[i]) From 4c950b831305002a7ab87e1d65bcb75b77eb143a Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Thu, 2 Jul 2026 12:46:28 +0200 Subject: [PATCH 18/19] Add fresh bm, python demo, and readme --- README.md | 1 + benchmark/CMakeLists.txt | 21 +++++ benchmark/bm_Fresh_L2Square.cpp | 153 ++++++++++++++++++++++++++++++++ demos/demo_Fresh_L2Square.py | 36 ++++++++ 4 files changed, 211 insertions(+) create mode 100644 benchmark/bm_Fresh_L2Square.cpp create mode 100644 demos/demo_Fresh_L2Square.py diff --git a/README.md b/README.md index eaf84b6..b7bce4c 100644 --- a/README.md +++ b/README.md @@ -47,6 +47,7 @@ The following table summarizes the key features of each algorithm: | **[SOFA](https://helios2.mi.parisdescartes.fr/~themisp/publications/icde25-sofa.pdf)** | In-memory similarity search using Symbolic Fourier Approximation (SFA) | | **[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) | diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt index a2eb474..f3e1d45 100644 --- a/benchmark/CMakeLists.txt +++ b/benchmark/CMakeLists.txt @@ -303,6 +303,27 @@ target_include_directories(bm_Hercules_L2Square ${CMAKE_CURRENT_SOURCE_DIR}/../commons ) +# ////// FRESH ////// +add_executable(bm_Fresh_L2Square + bm_Fresh_L2Square.cpp + bm_utils.cpp + ../commons/paramSetup.cpp + ../commons/test_bm_utils.cpp + ../commons/dataloaders.cpp +) +target_link_libraries(bm_Fresh_L2Square + PRIVATE + benchmark::benchmark + benchmark::benchmark_main + dino_lib + GTest::gtest +) +target_include_directories(bm_Fresh_L2Square + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../lib + ${CMAKE_CURRENT_SOURCE_DIR}/../commons +) + # ////// ODYSSEY (conditional on MPI) ////// if(BUILD_ODYSSEY_AVAILABLE) if(DEBUG_MSG) diff --git a/benchmark/bm_Fresh_L2Square.cpp b/benchmark/bm_Fresh_L2Square.cpp new file mode 100644 index 0000000..b7fdb62 --- /dev/null +++ b/benchmark/bm_Fresh_L2Square.cpp @@ -0,0 +1,153 @@ +#include +#include +#include +#include "bm_utils.hpp" +#include "../commons/dataloaders.hpp" +#include "../commons/VectorDataLoader.h" +#include "../commons/test_bm_utils.hpp" +#include "../lib/algos/Fresh.hpp" +#include "../lib/algos/DataSource.hpp" + +static bool endsWith(const std::string& s, const std::string& suffix) { + return s.size() >= suffix.size() && + s.compare(s.size() - suffix.size(), suffix.size(), suffix) == 0; +} + +struct FreshSearchOnlyFixture : public benchmark::Fixture { + daisy::Fresh* search = nullptr; + float* database = nullptr; + float* query = nullptr; + daisy::idx_t* I = nullptr; + float* D = nullptr; + daisy::idx_t n_query = 0; + size_t k = 0; + std::string dataset_name; + size_t n_database = 0; + int thread_count = 0; + + void SetUp(const benchmark::State& state) override { + int config_idx = static_cast(state.range(0)); + const SSTestConfig& config = test_configs_deep_seismic_astro270m[config_idx]; + + const bool use_fvecs = endsWith(config.dataset_path, ".fvecs") || endsWith(config.query_path, ".fvecs"); + size_t dim_u = 0, n_database_u = 0, n_q_u = 0; + database = nullptr; + + if (use_fvecs) { + database = fvecs_read(config.dataset_path.c_str(), &dim_u, &n_database_u, 0); + if (!database) { + std::cerr << "Failed to load dataset (fvecs)" << std::endl; + return; + } + const size_t query_limit = (config.query_limit > 0) ? static_cast(config.query_limit) : 0; + query = fvecs_read(config.query_path.c_str(), &dim_u, &n_q_u, query_limit); + if (!query) { + std::cerr << "Failed to load queries (fvecs)" << std::endl; + delete[] database; + return; + } + } else { + std::string dataset_filename = pathToFilename(config.dataset_path); + std::string query_filename = pathToFilename(config.query_path); + + daisy::idx_t dim, n_database, _, __; + if (!parseFilenameForConfig(dataset_filename, "bruteForce", dim, n_database, _, __)) { + std::cerr << "Failed to parse dataset config from filename: " << dataset_filename << std::endl; + return; + } + + daisy::idx_t dim_q, n_q, ___, ____; + if (!parseFilenameForConfig(query_filename, "bruteForce", dim_q, n_q, ___, ____)) { + std::cerr << "Failed to parse query config from filename: " << query_filename << std::endl; + return; + } + + if (dim != static_cast(dim_q)) { + std::cerr << "Dimension mismatch between dataset and queries" << std::endl; + return; + } + + dim_u = static_cast(dim); + n_database_u = static_cast(n_database); + if (config.query_limit > 0 && static_cast(config.query_limit) < n_q) + n_q = static_cast(config.query_limit); + n_q_u = static_cast(n_q); + + database = loadBinData(config.dataset_path.c_str(), n_database, dim, false); + if (!database) { + std::cerr << "Failed to load dataset" << std::endl; + return; + } + + query = loadBinData(config.query_path.c_str(), n_q, dim_q, false); + if (!query) { + std::cerr << "Failed to load queries" << std::endl; + delete[] database; + return; + } + } + + search = new daisy::Fresh(daisy::DistanceType::L2_SQUARED); + search->setNumThreads(config.thread_count); + search->setIndexWorkers(config.thread_count); + + fprintf(stderr, "[FRESH] Before buildIndex (n_database=%zu dim=%zu).\n", n_database_u, dim_u); + fflush(stderr); + + daisy::InMemoryDataSource data_source(database, static_cast(n_database_u), static_cast(dim_u)); + search->buildIndex(&data_source); + + fprintf(stderr, "[FRESH] Indexing finished (n_database=%zu dim=%zu).\n", n_database_u, dim_u); + fflush(stderr); + + k = static_cast(config.k_value); + n_query = static_cast(n_q_u); + I = new daisy::idx_t[n_query * k]; + D = new float[n_query * k]; + + dataset_name = config.name; + n_database = n_database_u; + thread_count = config.thread_count; + + fprintf(stderr, "[FRESH] n_database=%zu n_query=%zu dim=%zu k=%zu threads=%d\n", + n_database_u, (size_t)n_query, dim_u, k, config.thread_count); + fflush(stderr); + } + + void TearDown(const benchmark::State&) override { + fprintf(stderr, "[FRESH] TearDown start.\n"); + fflush(stderr); + delete search; + delete[] database; + delete[] query; + delete[] I; + delete[] D; + search = nullptr; + database = nullptr; + query = nullptr; + I = nullptr; + D = nullptr; + fprintf(stderr, "[FRESH] TearDown done.\n"); + fflush(stderr); + } +}; + +BENCHMARK_DEFINE_F(FreshSearchOnlyFixture, BM_Fresh_SearchOnly)(benchmark::State& state) { + for (auto _ : state) { + fprintf(stderr, "[FRESH] --- Query phase ---\n"); + fprintf(stderr, "[FRESH] dataset=%s n_database=%zu\n", dataset_name.c_str(), n_database); + fprintf(stderr, "[FRESH] search_threads=%d n_query=%zu k=%zu\n", thread_count, (size_t)n_query, k); + fflush(stderr); + search->searchIndex(query, n_query, static_cast(k), I, D); + fprintf(stderr, "[FRESH] Querying finished (n_query=%zu k=%zu).\n", (size_t)n_query, k); + fflush(stderr); + } +} + +BENCHMARK_REGISTER_F(FreshSearchOnlyFixture, BM_Fresh_SearchOnly) + // q=100, k=1,10,100,1000: DEEP (0-3), Seismic (4-7) + ->Args({0})->Args({1})->Args({2})->Args({3})->Args({4})->Args({5})->Args({6})->Args({7}) + ->Iterations(1) + ->Unit(benchmark::kMillisecond); + +BENCHMARK_MAIN(); diff --git a/demos/demo_Fresh_L2Square.py b/demos/demo_Fresh_L2Square.py new file mode 100644 index 0000000..6b39fb9 --- /dev/null +++ b/demos/demo_Fresh_L2Square.py @@ -0,0 +1,36 @@ +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, Fresh + +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 = Fresh(DistanceType.L2_SQUARED) + + index.setNumThreads(4) + index.buildIndex(db) + + 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() From b159bb8d90a2292778b3cfa8d3abb000ec63df1c Mon Sep 17 00:00:00 2001 From: ignaciomontoropineiro <138376065+ignaciomontoropineiro@users.noreply.github.com> Date: Thu, 2 Jul 2026 15:38:26 +0200 Subject: [PATCH 19/19] disable DumpyOS DTW tests --- tests/CMakeLists.txt | 52 ++++++++++++++++++++++---------------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 0c353c6..1202b97 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -653,32 +653,32 @@ gtest_discover_tests( WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} ) -# ////// DUMPYOS DTW ////// -add_executable( - test_DumpyOS_DTW - test_DumpyOS_DTW.cpp - test_utils.cpp -) - -target_link_libraries( - test_DumpyOS_DTW - PRIVATE - GTest::gtest_main - dino_lib - commons_lib - stdc++fs -) - -target_include_directories(test_DumpyOS_DTW - PRIVATE - ${CMAKE_CURRENT_SOURCE_DIR}/../lib - ${CMAKE_CURRENT_SOURCE_DIR}/../commons -) - -gtest_discover_tests( - test_DumpyOS_DTW - WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} -) +# # ////// DUMPYOS DTW ////// +# add_executable( +# test_DumpyOS_DTW +# test_DumpyOS_DTW.cpp +# test_utils.cpp +# ) +# +# target_link_libraries( +# test_DumpyOS_DTW +# PRIVATE +# GTest::gtest_main +# dino_lib +# commons_lib +# stdc++fs +# ) +# +# target_include_directories(test_DumpyOS_DTW +# PRIVATE +# ${CMAKE_CURRENT_SOURCE_DIR}/../lib +# ${CMAKE_CURRENT_SOURCE_DIR}/../commons +# ) +# +# gtest_discover_tests( +# test_DumpyOS_DTW +# WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} +# ) # ////// FRESH L2Square ////// add_executable(