Fix MESSI exact search with adaptive SAX and strided queries

README.md

@@ -78,6 +78,12 @@ SAX summary for each insert using the breakpoints established during the initial
 routes each new summary into the live iSAX tree and splits full leaves without rebuilding the
 index. Its first insert copies a borrowed initial in-memory database into owned growable storage.
 
+SAX encoding uses the ordered breakpoints for both tail intervals, including shifted and repeated
+adaptive cuts. `setNormalized(false)` selects adaptive breakpoints without changing input values.
+Python MESSI queries accept strided NumPy arrays (including `.fvecs` views); the binding copies
+queries to contiguous float32 storage when needed and checks the index dimension before searching.
+MESSI range refinement evaluates the full distance when a partial sum equals the inclusive radius.
+
 `FaissFlat` is the streaming baseline: `faiss::IndexFlatL2::add()` appends into flat storage with
 no training step and no rebuild, so nothing is frozen at build time that could drift from the live
 data and the index stays exact after every insert. It supports `L2_SQUARED` only.

docs/test-framework-guide.md

@@ -28,3 +28,43 @@ Follow the parameterized test pattern used throughout the test suite for BruteFo
 - Tests should be parameterized to run with different configurations (k-values, thread counts, datasets)
 - Include proper ground truth validation
 - Follow the existing naming conventions and code structure
+
+## MESSI exact top-k regressions
+
+`test_Messi_TopK_Distinct` compares every returned neighbor against exhaustive
+brute-force ground truth, checks distinct IDs, and recomputes the squared-L2
+distance of each returned vector. Missing neighbors are failures; a recall
+threshold below 100% is not sufficient for exact search. Integer-distance tie
+fixtures allow any distinct IDs with the exact expected distances.
+
+The parameterized cases cover Gaussian and adaptive breakpoints, one and four
+query workers, k=1/10/100/1000, repeated queries, and streaming updates. The
+inclusive-range regression checks zero and nonzero radii, exact boundary hits and
+partial sums equal to the radius before later nonzero chunks. The adaptive MESSI
+case in `test_equidepth` checks all k results with a small floating-point
+tolerance rather than only checking finite output. Simple CI
+runs these regressions through its complete CTest suite.
+
+`test_sax_breakpoints` covers positive/negative shifted cuts, repeated cuts, exact
+boundaries and both tail bins. It also checks scalar and SIMD record lower bounds
+against exact squared L2. The tests fail with the old sign-based SAX fallback.
+
+With `BUILD_PYTHON=ON`, `MessiPythonLayouts` additionally checks exact top-k and
+inclusive range results for contiguous, `.fvecs`-style, row-strided, Fortran and
+negative-stride queries, before and after streaming inserts. NumPy must be
+installed for this Python test. Python builds with Coconut disabled omit its binding.
+
+For the original unnormalized SIFT1M data, build Python bindings and install
+`numpy` and `faiss-cpu`, then run the optional manual audit:
+
+```bash
+PYTHONPATH=. python tests/check_messi_sift.py --dataset /path/to/sift --output sift-report.json
+```
+
+The audit uses 100 queries, k=1/10/100/1000, Gaussian/adaptive modes, one/four
+search workers, repeated queries, static 20K/1M builds and streaming 750K→875K→1M.
+It checks every distance against exact FAISS IndexFlatL2, recomputes returned
+ID/distance pairs, rejects duplicates and validates the full-dataset reference
+against the published SIFT ground truth. Distances of raw integer descriptors
+are compared exactly; equally distant alternative IDs are allowed. The dataset
+and FAISS are not required for the regular C++ tests.

lib/algos/Messi.cpp

@@ -128,7 +128,7 @@ namespace daisy
                 float dist = ts_euclidean_distance_SIMD(query, node->buffer->full_ts_buffer[i],
                                                         index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]);
 
-                if (dist <= pq_bsf->knn[pq_bsf->k - 1])
+                if (dist < pq_bsf->knn[pq_bsf->k - 1])
                 {
                     file_position_type pos = 0;
                     if (node->buffer->full_position_buffer != nullptr &&
@@ -137,7 +137,7 @@ namespace daisy
                         pos = *node->buffer->full_position_buffer[i] / index->settings->timeseries_size;
                     }
                     pthread_rwlock_wrlock(lock_queue);
-                    pqueue_bsf_insert(pq_bsf, dist, static_cast<long int>(pos), node);
+                    pqueue_bsf_insert_no_dup_position(pq_bsf, dist, static_cast<long int>(pos), node);
                     pthread_rwlock_unlock(lock_queue);
                 }
             }
@@ -147,7 +147,7 @@ namespace daisy
                 float dist = ts_euclidean_distance_SIMD(query, node->buffer->tmp_full_ts_buffer[i],
                                                         index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]);
 
-                if (dist <= pq_bsf->knn[pq_bsf->k - 1])
+                if (dist < pq_bsf->knn[pq_bsf->k - 1])
                 {
                     file_position_type pos = 0;
                     if (node->buffer->tmp_full_position_buffer != nullptr &&
@@ -156,7 +156,7 @@ namespace daisy
                         pos = *node->buffer->tmp_full_position_buffer[i] / index->settings->timeseries_size;
                     }
                     pthread_rwlock_wrlock(lock_queue);
-                    pqueue_bsf_insert(pq_bsf, dist, static_cast<long int>(pos), node);
+                    pqueue_bsf_insert_no_dup_position(pq_bsf, dist, static_cast<long int>(pos), node);
                     pthread_rwlock_unlock(lock_queue);
                 }
             }
@@ -174,11 +174,11 @@ namespace daisy
                     float dist = ts_euclidean_distance_SIMD(query, &(rawfile[*node->buffer->partial_position_buffer[i]]),
                                                             index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]);
 
-                    if (dist <= pq_bsf->knn[pq_bsf->k - 1])
+                    if (dist < pq_bsf->knn[pq_bsf->k - 1])
                     {
                         pthread_rwlock_wrlock(lock_queue);
 
-                        pqueue_bsf_insert(pq_bsf, dist, *node->buffer->partial_position_buffer[i] / index->settings->timeseries_size, node);
+                        pqueue_bsf_insert_no_dup_position(pq_bsf, dist, *node->buffer->partial_position_buffer[i] / index->settings->timeseries_size, node);
 
                         pthread_rwlock_unlock(lock_queue);
                     }
@@ -229,7 +229,7 @@ namespace daisy
                 if (dist <= pq_bsf->knn[pq_bsf->k - 1])
                 {
                     pthread_rwlock_wrlock(lock_queue);
-                    pqueue_bsf_insert(pq_bsf, dist, *node->buffer->partial_position_buffer[i] / index->settings->timeseries_size, node);
+                    pqueue_bsf_insert_no_dup_position(pq_bsf, dist, *node->buffer->partial_position_buffer[i] / index->settings->timeseries_size, node);
                     pthread_rwlock_unlock(lock_queue);
                 }
             }
@@ -496,10 +496,13 @@ namespace daisy
         if (node == NULL || node->buffer == NULL)
             return;
 
+        // The bounded SIMD kernel stops at >= bound. Inclusive ranges must
+        // continue past partial sums equal to r before accepting a full distance.
+        const float evaluation_bound = std::nextafter(r, std::numeric_limits<float>::infinity());
         for (int i = 0; i < node->buffer->full_buffer_size; i++)
         {
             float dist = ts_euclidean_distance_SIMD(query, node->buffer->full_ts_buffer[i],
-                                                    index->settings->timeseries_size, r);
+                                                    index->settings->timeseries_size, evaluation_bound);
             if (dist <= r)
             {
                 file_position_type pos = 0;
@@ -515,7 +518,7 @@ namespace daisy
         for (int i = 0; i < node->buffer->tmp_full_buffer_size; i++)
         {
             float dist = ts_euclidean_distance_SIMD(query, node->buffer->tmp_full_ts_buffer[i],
-                                                    index->settings->timeseries_size, r);
+                                                    index->settings->timeseries_size, evaluation_bound);
             if (dist <= r)
             {
                 file_position_type pos = 0;
@@ -540,7 +543,7 @@ namespace daisy
             if (distmin <= r)
             {
                 float dist = ts_euclidean_distance_SIMD(query, &(rawfile[*node->buffer->partial_position_buffer[i]]),
-                                                        index->settings->timeseries_size, r);
+                                                        index->settings->timeseries_size, evaluation_bound);
                 if (dist <= r)
                 {
                     pthread_rwlock_wrlock(lock);
@@ -558,6 +561,8 @@ namespace daisy
         ts_type *paa = wd->paa;
         ts_type *ts = wd->ts;
         float r = wd->r;
+        // The shared top-k enqueue helper uses <; range pruning is inclusive.
+        const float pruning_bound = std::nextafter(r, std::numeric_limits<float>::infinity());
         int n_pqueue = wd->n_pqueue;
         float *rawfile = wd->rawfile;
         int tnumber = rand() % n_pqueue;
@@ -570,7 +575,7 @@ namespace daisy
             int current_root_node_number = __sync_fetch_and_add(wd->node_counter, 1);
             if (current_root_node_number >= wd->amountnode)
                 break;
-            insert_tree_node_m_hybridpqueue(paa, wd->nodelist[current_root_node_number], index, r,
+            insert_tree_node_m_hybridpqueue(paa, wd->nodelist[current_root_node_number], index, pruning_bound,
                                             wd->allpq, wd->alllock, &tnumber, n_pqueue);
         }
 

lib/isax/SAX.cpp

@@ -1,5 +1,6 @@
 #include <math.h>
 #include <float.h>
+#include <algorithm>
 
 #include "SAX.hpp"
 #include "../utils/SIMD.hpp"
@@ -86,22 +87,13 @@ namespace daisy
         //       TO   (c - 1) * (c - 2) / 2 + c - 1
         int offset = ((cardinality - 1) * (cardinality - 2)) / 2;
 
-        int si;
-        for (si = 0; si < segments; si++)
-        {
-            sax_out[si] = 0;
-
-            // First object = daisy_active_breakpoints[offset]
-            // Last object = daisy_active_breakpoints[offset + cardinality - 2]
-            // Size of sub-array = cardinality - 1
-
-            float *res = (float *)bsearch(&paa[si], &daisy_active_breakpoints[offset], cardinality - 1,
-                                          sizeof(ts_type), compare);
-            if (res != NULL)
-                sax_out[si] = (int)(res - &daisy_active_breakpoints[offset]);
-            else if (paa[si] > 0)
-                sax_out[si] = cardinality - 1;
-        }
+        // Each interval is (previous_cut, cut], with unbounded tail bins.
+        // Searching the cuts directly also handles shifted/repeated adaptive
+        // cuts; neither the preceding cardinality nor the sign is an endpoint.
+        const float *cuts = daisy_active_breakpoints + offset;
+        for (int si = 0; si < segments; si++)
+            sax_out[si] = static_cast<sax_type>(
+                std::lower_bound(cuts, cuts + cardinality - 1, paa[si]) - cuts);
 
         free(paa);
         return SUCCESS;
@@ -126,22 +118,10 @@ namespace daisy
     {
         int offset = ((cardinality - 1) * (cardinality - 2)) / 2;
 
-        int si;
-        for (si = 0; si < segments; si++)
-        {
-            sax[si] = 0;
-
-            // First object = daisy_active_breakpoints[offset]
-            // Last object = daisy_active_breakpoints[offset + cardinality - 2]
-            // Size of sub-array = cardinality - 1
-
-            float *res = (float *)bsearch(&paa[si], &daisy_active_breakpoints[offset], cardinality - 1,
-                                          sizeof(ts_type), compare);
-            if (res != NULL)
-                sax[si] = (int)(res - &daisy_active_breakpoints[offset]);
-            else if (paa[si] > 0)
-                sax[si] = cardinality - 1;
-        }
+        const float *cuts = daisy_active_breakpoints + offset;
+        for (int si = 0; si < segments; si++)
+            sax[si] = static_cast<sax_type>(
+                std::lower_bound(cuts, cuts + cardinality - 1, paa[si]) - cuts);
 
         return SUCCESS;
     }

lib/isax/iSAXPqueue.cpp

@@ -498,11 +498,26 @@ namespace daisy
     {
         int i, j;
 
-        // Fast duplicate check: if position already present, skip
+        // If the position is already queued, keep only its best distance. An
+        // earlier entry may carry an early-abandoned (over-estimated) value, so
+        // a better distance replaces it instead of being dropped.
         for (i = 0; i < q->k; i++)
         {
             if (q->position[i] == position)
-                return;
+            {
+                if (data >= q->knn[i])
+                    return;
+                for (j = i; j < q->k - 1; j++)
+                {
+                    q->knn[j] = q->knn[j + 1];
+                    q->position[j] = q->position[j + 1];
+                    q->node[j] = q->node[j + 1];
+                }
+                q->knn[q->k - 1] = FLT_MAX;
+                q->position[q->k - 1] = -1;
+                q->node[q->k - 1] = nullptr;
+                break;
+            }
         }
 
         for (i = 0; i < q->k; i++)

lib/isax/iSAXSearch.cpp

@@ -22,7 +22,7 @@ namespace daisy
             {
                 float dist = ts_euclidean_distance_SIMD(query, node->buffer->full_ts_buffer[i],
                                                         index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]);
-                if (dist <= pq_bsf->knn[pq_bsf->k - 1])
+                if (dist < pq_bsf->knn[pq_bsf->k - 1])
                 {
                     file_position_type pos = 0;
                     if (node->buffer->full_position_buffer != nullptr &&
@@ -36,7 +36,7 @@ namespace daisy
             {
                 float dist = ts_euclidean_distance_SIMD(query, node->buffer->tmp_full_ts_buffer[i],
                                                         index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]);
-                if (dist <= pq_bsf->knn[pq_bsf->k - 1])
+                if (dist < pq_bsf->knn[pq_bsf->k - 1])
                 {
                     file_position_type pos = 0;
                     if (node->buffer->tmp_full_position_buffer != nullptr &&
@@ -53,7 +53,7 @@ namespace daisy
             float dist = ts_euclidean_distance_SIMD(query, &(rawfile[*node->buffer->partial_position_buffer[i]]),
                                                     index->settings->timeseries_size, pq_bsf->knn[pq_bsf->k - 1]);
 
-            if (dist <= pq_bsf->knn[pq_bsf->k - 1])
+            if (dist < pq_bsf->knn[pq_bsf->k - 1])
             {
                 file_position_type pos = *node->buffer->partial_position_buffer[i] / index->settings->timeseries_size;
                 pqueue_bsf_insert_no_dup_position(pq_bsf, dist, pos, node);

pybinds/setup.cpp

@@ -12,7 +12,9 @@
 #include "../lib/distance_computers/DistanceComputer.hpp"
 #include "../lib/algos/Bruteforce.hpp"
 #include "../lib/algos/LbBruteforce.hpp"
+#ifdef COCONUT_ENABLED
 #include "../lib/algos/Coconut.hpp"
+#endif
 #include "../lib/algos/Messi.hpp"
 #if ODYSSEY_MPI
 #include "../lib/algos/hodyssey/Odyssey.hpp"
@@ -377,6 +379,7 @@ PYBIND11_MODULE(_core, m)
             return pybind11::make_tuple(I, D); }, "Search using SearchConfig (top-k or range) and return (indices, distances)");
 
     ////// COCONUT //////
+#ifdef COCONUT_ENABLED
     pybind11::class_<daisy::Coconut>(m, "Coconut", "COCONUT sortable-SAX index (static bottom-up build + streaming insert)")
         .def(pybind11::init<daisy::DistanceType>(), "Create a new Coconut with the given distance metric")
 
@@ -447,6 +450,8 @@ PYBIND11_MODULE(_core, m)
             self.searchIndex(static_cast<float *>(query_buf.ptr), n_query, config, I, D);
             return pybind11::make_tuple(I, D); }, "Search using SearchConfig (top-k or range) and return (indices, distances)");
 
+#endif
+
     ////// MESSI //////
     pybind11::class_<daisy::Messi>(m, "Messi", "MESSI (Multi-Queue Efficient SAX Similarity Index) algorithm for time series similarity search")
         // Constructor
@@ -534,16 +539,22 @@ PYBIND11_MODULE(_core, m)
              "Incrementally insert a batch of series into the live MESSI index")
 
         // Search the index with query array and return top-k results
-        .def("searchIndex", [](daisy::Messi &self, pybind11::array_t<float> query, daisy::idx_t k)
+        .def("searchIndex", [](daisy::Messi &self,
+                               pybind11::array_t<float, pybind11::array::c_style |
+                                                       pybind11::array::forcecast> query,
+                               daisy::idx_t k)
              {
             pybind11::buffer_info query_buf = query.request();
             if (query_buf.ndim != 2)
                 throw std::runtime_error("Query array must be 2D");
+            if (self.getDim() == 0 || self.getNDatabase() == 0)
+                throw std::runtime_error("MESSI index must be built before searching");
+            if (static_cast<daisy::idx_t>(query_buf.shape[1]) != self.getDim())
+                throw std::runtime_error("Query dimension does not match the index dimension");
             if (k <= 0)
                 throw std::runtime_error("k must be positive");
 
             const daisy::idx_t n_query = query_buf.shape[0];
-            const daisy::idx_t dim = query_buf.shape[1];
 
             std::vector<daisy::idx_t> indices(n_query * k);
             std::vector<float> distances(n_query * k);
@@ -554,11 +565,20 @@ PYBIND11_MODULE(_core, m)
                 pybind11::array_t<daisy::idx_t>({n_query, k}, indices.data()),
                 pybind11::array_t<float>({n_query, k}, distances.data())
             ); }, "Search the MESSI index using queries and return (indices, distances)")
-        .def("searchIndex", [](daisy::Messi &self, pybind11::array_t<float> query, daisy::SearchConfig config)
+        .def("searchIndex", [](daisy::Messi &self,
+                               pybind11::array_t<float, pybind11::array::c_style |
+                                                       pybind11::array::forcecast> query,
+                               daisy::SearchConfig config)
              {
             pybind11::buffer_info query_buf = query.request();
             if (query_buf.ndim != 2)
                 throw std::runtime_error("Query array must be 2D");
+            if (self.getDim() == 0 || self.getNDatabase() == 0)
+                throw std::runtime_error("MESSI index must be built before searching");
+            if (static_cast<daisy::idx_t>(query_buf.shape[1]) != self.getDim())
+                throw std::runtime_error("Query dimension does not match the index dimension");
+            if (config.type == daisy::QueryType::TOP_K && config.k == 0)
+                throw std::runtime_error("k must be positive");
             daisy::idx_t n_query = query_buf.shape[0];
             std::vector<std::vector<daisy::idx_t>> I;
             std::vector<std::vector<float>> D;