From 1ce96bb8ac7f5ec8ad78cf7952c7fc2c6f03bc41 Mon Sep 17 00:00:00 2001 From: Kristofferjoha Date: Wed, 28 Oct 2026 12:43:44 +0100 Subject: [PATCH 1/4] heap and lower check improvements --- engine/src/query/knn.cpp | 37 ++++++++++++++++++++++++++----------- engine/tests/knn_test.cpp | 13 +++++++++++++ 2 files changed, 39 insertions(+), 11 deletions(-) diff --git a/engine/src/query/knn.cpp b/engine/src/query/knn.cpp index 90d9339..3b4a21a 100644 --- a/engine/src/query/knn.cpp +++ b/engine/src/query/knn.cpp @@ -192,25 +192,38 @@ std::vector knn(const core::TrajectoryDb &trajectory_db, throw std::invalid_argument("knn: k = 0, no neighbors"); } - // Lower bound on the Fréchet distance from the query to each trajectory, sorted ascending. + // Extract candidates incrementally by the stronger of the endpoint and box-extent bounds. // Once a bound reaches the current k-th best distance, every remaining candidate is out too. std::vector candidates(trajectory_db.size()); for (std::size_t i = 0; i < trajectory_db.size(); ++i) { const core::TrajectoryView reference = trajectory_db.trajectory(i); - candidates[i] = {squared_endpoint_lower_bound(query_trajectory, reference), i}; + candidates[i] = { + std::max(squared_endpoint_lower_bound(query_trajectory, reference), + squared_box_extent_lower_bound(query_trajectory, reference)), + i}; } - std::ranges::sort(candidates, [](const Candidate &a, const Candidate &b) { - return std::tie(a.squared_lower_bound, a.idx) < std::tie(b.squared_lower_bound, b.idx); - }); + const auto farther = [](const Candidate &a, const Candidate &b) { + return std::tie(a.squared_lower_bound, a.idx) > std::tie(b.squared_lower_bound, b.idx); + }; + std::make_heap(candidates.begin(), candidates.end(), farther); + const auto extract_candidate = [&]() { + std::pop_heap(candidates.begin(), candidates.end(), farther); + const Candidate candidate = candidates.back(); + candidates.pop_back(); + return candidate; + }; double initial_cutoff = std::numeric_limits::infinity(); + std::vector seed_candidates; if (k <= candidates.size()) { constexpr std::size_t kInitialBoundCandidates = 64; const auto count = std::min(candidates.size(), std::max(k, kInitialBoundCandidates)); + seed_candidates.reserve(count); std::priority_queue upper_bounds; for (std::size_t i = 0; i < count; ++i) { + seed_candidates.push_back(extract_candidate()); const double upper = squared_greedy_upper_bound( - query_trajectory, trajectory_db.trajectory(candidates[i].idx)); + query_trajectory, trajectory_db.trajectory(seed_candidates.back().idx)); if (upper_bounds.size() < k) upper_bounds.push(upper); else if (upper < upper_bounds.top()) { @@ -226,7 +239,11 @@ std::vector knn(const core::TrajectoryDb &trajectory_db, std::priority_queue curr_neighbors; std::vector prev; std::vector curr; - for (const Candidate &candidate : candidates) { + std::size_t seed_index = 0; + while (seed_index < seed_candidates.size() || !candidates.empty()) { + const Candidate candidate = seed_index < seed_candidates.size() + ? seed_candidates[seed_index++] + : extract_candidate(); const std::size_t i = candidate.idx; // The greedy paths bound the k-th distance until exact results tighten it. @@ -241,10 +258,8 @@ std::vector knn(const core::TrajectoryDb &trajectory_db, core::TrajectoryView reference_trajectory = trajectory_db.trajectory(i); if (std::isfinite(curr_max_dist) && - (squared_box_extent_lower_bound(query_trajectory, reference_trajectory) >= - curr_max_dist || - squared_points_to_box_lower_bound(query_trajectory, reference_trajectory, - curr_max_dist) >= curr_max_dist)) { + squared_points_to_box_lower_bound(query_trajectory, reference_trajectory, + curr_max_dist) >= curr_max_dist) { continue; } double curr_dist = get_squared_frechet_distance(query_trajectory, reference_trajectory, diff --git a/engine/tests/knn_test.cpp b/engine/tests/knn_test.cpp index 4de26d5..9f9a856 100644 --- a/engine/tests/knn_test.cpp +++ b/engine/tests/knn_test.cpp @@ -148,6 +148,19 @@ TEST(Knn, MatchesBruteForceForOverlappingLongTrajectories) { } } +TEST(Knn, MatchesBruteForceAcrossSeedCountBoundary) { + std::mt19937_64 rng(31); + const auto db = make_db(random_tracks(rng, 90)); + const auto queries = make_db(random_tracks(rng, 3)); + for (const std::size_t k : {1u, 63u, 64u, 65u, 90u, 100u}) { + for (std::size_t q = 0; q < queries.size(); ++q) { + EXPECT_EQ(query::knn(db, queries.trajectory(q), k), + brute_force_knn(db, queries.trajectory(q), k)) + << "query " << q << ", k " << k; + } + } +} + TEST(Knn, AcceptsQueryViewWithoutCachedBoundingBox) { std::mt19937_64 rng(23); const auto db = make_db(random_tracks(rng, 30)); From 406bcb1745bc81f36995d97777c0ec7b4b2ab1f9 Mon Sep 17 00:00:00 2001 From: TB-PC Date: Thu, 29 Oct 2026 00:07:04 +0100 Subject: [PATCH 2/4] init --- engine/src/query/knn.cpp | 44 +++++++++++++++++++++++++++++----------- 1 file changed, 32 insertions(+), 12 deletions(-) diff --git a/engine/src/query/knn.cpp b/engine/src/query/knn.cpp index 3b4a21a..069e7e7 100644 --- a/engine/src/query/knn.cpp +++ b/engine/src/query/knn.cpp @@ -52,28 +52,48 @@ double get_squared_frechet_distance(const core::TrajectoryView &query, if (last == 0) return inf; + // Cells at or above the cutoff keep their value instead of becoming infinity. They can only + // lead to values at or above the cutoff, so every value below it stays exact. std::size_t first = 0; for (std::size_t i = 1; i < n; ++i) { + const double qx = query.x[i]; + const double qy = query.y[i]; + const auto dist = [&](std::size_t j) { + return squared_point_dist(qx, qy, reference.x[j], reference.y[j]); + }; std::size_t next_first = m; std::size_t next_last = 0; - double left = inf; - for (std::size_t j = first; j < m; ++j) { - // Beyond the previous row, only a reachable cell on the left can - // extend the path. j == last can still use the diagonal predecessor. - if (j > last && left >= max_squared_dist) - break; - const double above = j < last ? prev[j] : inf; - const double diagonal = j > first && j <= last ? prev[j - 1] : inf; - const double predecessor = std::min({above, diagonal, left}); - value = std::max(predecessor, squared_point_dist(query.x[i], query.y[i], reference.x[j], - reference.y[j])); - left = value < max_squared_dist ? value : inf; + + // Within the previous row's reachable span [first, last), the cell above always exists, + // and the diagonal does for every cell but the first + double left = std::max(prev[first], dist(first)); + curr[first] = left; + if (left < max_squared_dist) { + next_first = first; + next_last = first + 1; + } + for (std::size_t j = first + 1; j < last; ++j) { + left = std::max(std::min({prev[j], prev[j - 1], left}), dist(j)); curr[j] = left; if (left < max_squared_dist) { next_first = std::min(next_first, j); next_last = j + 1; } } + + // Past the span only the cell on the left continues a path, plus the diagonal at j == last + for (std::size_t j = last; j < m; ++j) { + const double predecessor = j == last ? std::min(prev[last - 1], left) : left; + if (predecessor >= max_squared_dist) + break; + left = std::max(predecessor, dist(j)); + if (left >= max_squared_dist) + break; + curr[j] = left; + next_first = std::min(next_first, j); + next_last = j + 1; + } + if (next_last == 0) return inf; first = next_first; From 46bb1a84747967c524f0fa182483ab47c23aa4ff Mon Sep 17 00:00:00 2001 From: Kristofferjoha Date: Thu, 29 Oct 2026 13:09:26 +0100 Subject: [PATCH 3/4] format --- engine/src/query/knn.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/engine/src/query/knn.cpp b/engine/src/query/knn.cpp index 069e7e7..790fc2d 100644 --- a/engine/src/query/knn.cpp +++ b/engine/src/query/knn.cpp @@ -217,10 +217,9 @@ std::vector knn(const core::TrajectoryDb &trajectory_db, std::vector candidates(trajectory_db.size()); for (std::size_t i = 0; i < trajectory_db.size(); ++i) { const core::TrajectoryView reference = trajectory_db.trajectory(i); - candidates[i] = { - std::max(squared_endpoint_lower_bound(query_trajectory, reference), - squared_box_extent_lower_bound(query_trajectory, reference)), - i}; + candidates[i] = {std::max(squared_endpoint_lower_bound(query_trajectory, reference), + squared_box_extent_lower_bound(query_trajectory, reference)), + i}; } const auto farther = [](const Candidate &a, const Candidate &b) { return std::tie(a.squared_lower_bound, a.idx) > std::tie(b.squared_lower_bound, b.idx); From ce11e74f569aeca7d5b3e99b8dc4ca580c9b11d7 Mon Sep 17 00:00:00 2001 From: Kristofferjoha Date: Thu, 29 Oct 2026 13:49:45 +0100 Subject: [PATCH 4/4] two minor speedups --- engine/src/query/where_at_query.cpp | 12 ++++++++++++ engine/src/selection/rdp.cpp | 11 ++++++----- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/engine/src/query/where_at_query.cpp b/engine/src/query/where_at_query.cpp index 3031880..fcfabee 100644 --- a/engine/src/query/where_at_query.cpp +++ b/engine/src/query/where_at_query.cpp @@ -25,6 +25,7 @@ void validate_config(const WhereAtConfig &config) { std::vector eligible_intervals(const core::TrajectoryView &trajectory, double max_gap_seconds) { std::vector intervals; + intervals.reserve(trajectory.size() > 0 ? trajectory.size() - 1 : 0); double duration = 0; for (std::size_t j = 1; j < trajectory.size(); ++j) { const double start = static_cast(trajectory.t[j - 1]); @@ -65,6 +66,17 @@ WhereAtWorkload generate_where_at_queries(const core::TrajectoryDb &uncompressed validate_config(config); WhereAtWorkload result; + std::size_t eligible = 0; + for (std::size_t i = 0; i < uncompressed_database.size(); ++i) { + const auto times = uncompressed_database.trajectory(i).t; + eligible += + std::adjacent_find(times.begin(), times.end(), [&](auto a, auto b) { + return static_cast(b) - static_cast(a) <= config.max_gap_seconds; + }) != times.end(); + } + if (eligible > result.queries.max_size() / config.probes_per_trajectory) + throw std::length_error("where-at: too many probes"); + result.queries.reserve(eligible * config.probes_per_trajectory); std::mt19937_64 rng(seed); for (std::size_t i = 0; i < uncompressed_database.size(); i++) { const auto trajectory = uncompressed_database.trajectory(i); diff --git a/engine/src/selection/rdp.cpp b/engine/src/selection/rdp.cpp index 4a041ae..1543612 100644 --- a/engine/src/selection/rdp.cpp +++ b/engine/src/selection/rdp.cpp @@ -1,7 +1,6 @@ #include "backend/engine/selection/rdp.hpp" #include -#include #include #include @@ -14,14 +13,16 @@ namespace { // The implementation uses max-prio queue with fixed point budget instead of standard epsilon error // like in vid. -double segment_distance(double px, double py, double ax, double ay, double bx, double by) { +double squared_segment_distance(double px, double py, double ax, double ay, double bx, double by) { const double dx = bx - ax; const double dy = by - ay; const double len_sq = dx * dx + dy * dy; double s = 0.0; if (len_sq > 0.0) s = std::clamp(((px - ax) * dx + (py - ay) * dy) / len_sq, 0.0, 1.0); - return std::hypot(px - (ax + s * dx), py - (ay + s * dy)); + const double ex = px - (ax + s * dx); + const double ey = py - (ay + s * dy); + return ex * ex + ey * ey; } struct Split { @@ -39,7 +40,7 @@ void push_split(std::priority_queue &heap, std::span x, return; Split split{-1.0, first + 1, first, last}; for (std::size_t k = first + 1; k < last; ++k) { - const double d = segment_distance(x[k], y[k], x[first], y[first], x[last], y[last]); + const double d = squared_segment_distance(x[k], y[k], x[first], y[first], x[last], y[last]); if (d > split.error) { split.error = d; split.index = k; @@ -87,4 +88,4 @@ std::vector simplify(const core::TrajectoryView &trajectory, return indices; } -} // namespace selection::rdp \ No newline at end of file +} // namespace selection::rdp