diff --git a/engine/configs/geolife.json b/engine/configs/geolife.json index e550a0f..a6aa57b 100644 --- a/engine/configs/geolife.json +++ b/engine/configs/geolife.json @@ -1,5 +1,5 @@ { - "data_path": "output/geolife.parquet", + "data_path": "data/geolife.parquet", "seed": 42, "budget_shares": [0.0025, 0.005, 0.01, 0.02], "split": { @@ -8,8 +8,8 @@ }, "algorithms": ["RDP", "SED"], "knn": { - "k": 5, - "query_count": 10 + "k": 10, + "query_count": 0 }, "range": { "query_count": 500, diff --git a/engine/include/backend/engine/core/trajectory_db.hpp b/engine/include/backend/engine/core/trajectory_db.hpp index 8ed1c3f..c35cff6 100644 --- a/engine/include/backend/engine/core/trajectory_db.hpp +++ b/engine/include/backend/engine/core/trajectory_db.hpp @@ -3,11 +3,20 @@ #include #include #include +#include #include #include namespace core { +// Smallest axis-aligned box around a trajectory's x/y. +struct BoundingBox { + double x_min = -std::numeric_limits::infinity(); + double x_max = std::numeric_limits::infinity(); + double y_min = -std::numeric_limits::infinity(); + double y_max = std::numeric_limits::infinity(); +}; + struct TrajectoryView { std::int32_t trajectory_id; std::span t; @@ -15,6 +24,7 @@ struct TrajectoryView { std::span lon; std::span x; std::span y; + BoundingBox bbox{}; std::size_t size() const { return t.size(); } }; @@ -32,8 +42,13 @@ class TrajectoryDb { const auto start = m_offsets[i]; const auto len = m_offsets[i + 1] - start; return TrajectoryView{ - m_trajectory_ids[i], {m_t.data() + start, len}, {m_lat.data() + start, len}, - {m_lon.data() + start, len}, {m_x.data() + start, len}, {m_y.data() + start, len}, + m_trajectory_ids[i], + {m_t.data() + start, len}, + {m_lat.data() + start, len}, + {m_lon.data() + start, len}, + {m_x.data() + start, len}, + {m_y.data() + start, len}, + m_boxes[i], }; } @@ -48,6 +63,7 @@ class TrajectoryDb { std::vector m_lon; std::vector m_x; std::vector m_y; + std::vector m_boxes; }; // Copies validated trajectories, including projected coordinates, in the given order. diff --git a/engine/include/backend/engine/evaluator/knn.hpp b/engine/include/backend/engine/evaluator/knn.hpp index 981f38a..2c3664b 100644 --- a/engine/include/backend/engine/evaluator/knn.hpp +++ b/engine/include/backend/engine/evaluator/knn.hpp @@ -12,6 +12,8 @@ namespace evaluator { struct KnnScores { double kendalls_tau = std::numeric_limits::quiet_NaN(); double recall_at_k = std::numeric_limits::quiet_NaN(); + std::size_t query_count = 0; + std::size_t tau_query_count = 0; // queries with at least two shared neighbors }; // Each score is averaged over the queries where it is defined, NaN if it is defined for none. @@ -21,6 +23,7 @@ KnnScores evaluate_knn(const core::TrajectoryDb &compressed, const query::KnnWor namespace evaluator::knn { +// Rank agreement on shared neighbors only; NaN when fewer than two are shared. double kendalls_tau(const std::vector &original, const std::vector &simplified); diff --git a/engine/include/backend/engine/query/knn.hpp b/engine/include/backend/engine/query/knn.hpp index ed9e7ed..9116d9b 100644 --- a/engine/include/backend/engine/query/knn.hpp +++ b/engine/include/backend/engine/query/knn.hpp @@ -31,7 +31,7 @@ KnnWorkload generate_knn_queries(const core::TrajectoryDb &original, std::uint64_t seed = 42); // Indices into trajectories of the k trajectories closest to query by discrete Fréchet -// distance, closest first. +// distance, closest first. Equal distances are ordered by database index. std::vector knn(const core::TrajectoryDb &trajectories, const core::TrajectoryView &query, std::size_t k); diff --git a/engine/src/core/trajectory_db.cpp b/engine/src/core/trajectory_db.cpp index 05da944..3fd5ce3 100644 --- a/engine/src/core/trajectory_db.cpp +++ b/engine/src/core/trajectory_db.cpp @@ -1,5 +1,7 @@ #include "backend/engine/core/trajectory_db.hpp" #include "backend/engine/core/projection.hpp" +#include +#include #include #include #include @@ -37,6 +39,21 @@ TrajectoryDb::TrajectoryDb(std::vector offsets, xs[k] = p.x; ys[k] = p.y; } + + // For each trajectory it goes throguh all its points + // and keeps min/max of x and y coordinates seen + constexpr double inf = std::numeric_limits::infinity(); + m_boxes.resize(size()); + for (std::size_t i = 0; i < size(); ++i) { + BoundingBox box{inf, -inf, inf, -inf}; + for (auto k = m_offsets[i]; k < m_offsets[i + 1]; ++k) { + box.x_min = std::min(box.x_min, xs[k]); + box.x_max = std::max(box.x_max, xs[k]); + box.y_min = std::min(box.y_min, ys[k]); + box.y_max = std::max(box.y_max, ys[k]); + } + m_boxes[i] = box; + } } TrajectoryDb subset(const TrajectoryDb &trajectories, std::span indices) { @@ -55,10 +72,12 @@ TrajectoryDb subset(const TrajectoryDb &trajectories, std::span(tau_scored); if (recall_scored != 0) diff --git a/engine/src/pipeline/pipeline.cpp b/engine/src/pipeline/pipeline.cpp index f96a79c..a50befd 100644 --- a/engine/src/pipeline/pipeline.cpp +++ b/engine/src/pipeline/pipeline.cpp @@ -45,7 +45,8 @@ std::size_t total_budget(const core::TrajectoryDb &trajectories, double share) { } void print_scores(const evaluator::Scoretable &scores) { - std::cout << " knn tau: " << scores.knn.kendalls_tau + std::cout << " knn tau on shared neighbors: " << scores.knn.kendalls_tau << " (" + << scores.knn.tau_query_count << '/' << scores.knn.query_count << " queries)" << ", recall@k: " << scores.knn.recall_at_k << '\n'; std::cout << " where-at: within radius=" << scores.where_at.within_radius << ", average hit=" << scores.where_at.average_hit << '\n'; @@ -57,10 +58,11 @@ void print_scores(const evaluator::Scoretable &scores) { void run(const Config &config) { const auto [database, queries] = load_split(config); print_size("database", database); - print_size("queries", queries); + print_size("query pool", queries); const auto knn_workload = query::generate_knn_queries(database, queries, config.knn, config.seed); + std::cout << "knn: " << knn_workload.entries.size() << " queries, k=" << knn_workload.k << '\n'; const auto where_at = query::generate_where_at_queries(database, config.where_at, config.seed); std::cout << "where-at: " << where_at.queries.size() << " probes, " << where_at.skipped_trajectories << " trips without eligible time\n"; diff --git a/engine/src/query/knn.cpp b/engine/src/query/knn.cpp index 804087f..6c8d414 100644 --- a/engine/src/query/knn.cpp +++ b/engine/src/query/knn.cpp @@ -2,12 +2,14 @@ #include "backend/engine/core/parallel.hpp" #include "backend/engine/core/trajectory_db.hpp" #include +#include #include #include #include #include #include #include +#include #include namespace query { @@ -33,29 +35,52 @@ double get_squared_frechet_distance(const core::TrajectoryView &query, prev.resize(m); curr.resize(m); - for (std::size_t i = 0; i < n; ++i) { - double row_min = std::numeric_limits::infinity(); - for (std::size_t j = 0; j < m; ++j) { - const double d = - squared_point_dist(query.x[i], query.y[i], reference.x[j], reference.y[j]); - if (i == 0 && j == 0) { - curr[j] = d; - } else if (i == 0) { - curr[j] = std::max(d, curr[j - 1]); - } else if (j == 0) { - curr[j] = std::max(d, prev[j]); - } else { - curr[j] = std::max(d, std::min({prev[j], curr[j - 1], prev[j - 1]})); + const double inf = std::numeric_limits::infinity(); + + // Only cells below the current k-th distance can improve the answer. Track the + // reachable span in each row instead of evaluating the entire n-by-m grid. + double value = 0.0; + std::size_t last = 0; + for (std::size_t j = 0; j < m; ++j) { + value = std::max( + value, squared_point_dist(query.x[0], query.y[0], reference.x[j], reference.y[j])); + if (value >= max_squared_dist) + break; + prev[j] = value; + last = j + 1; + } + if (last == 0) + return inf; + + std::size_t first = 0; + for (std::size_t i = 1; i < n; ++i) { + 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; + curr[j] = left; + if (left < max_squared_dist) { + next_first = std::min(next_first, j); + next_last = j + 1; } - row_min = std::min(row_min, curr[j]); - } - // Every path to the end passes through this row, so the result is at least row_min. - if (row_min >= max_squared_dist) { - return row_min; } + if (next_last == 0) + return inf; + first = next_first; + last = next_last; std::swap(prev, curr); } - return prev[m - 1]; + return last == m ? prev[m - 1] : inf; } double squared_endpoint_lower_bound(const core::TrajectoryView &query, @@ -69,6 +94,86 @@ double squared_endpoint_lower_bound(const core::TrajectoryView &query, squared_point_dist(query.x[q_last], query.y[q_last], reference.x[r_last], reference.y[r_last])); } + +// Point distance to bounding rectangle of a given trajectory bounding box +double squared_dist_to_box(double x, double y, const core::BoundingBox &box) { + const double dx = std::max({box.x_min - x, 0.0, x - box.x_max}); + const double dy = std::max({box.y_min - y, 0.0, y - box.y_max}); + return dx * dx + dy * dy; +} + +double squared_box_extent_lower_bound(const core::TrajectoryView &query, + const core::TrajectoryView &reference) { + if (query.size() == 0 || reference.size() == 0) { + return std::numeric_limits::infinity(); + } + const core::BoundingBox &a = query.bbox; + const core::BoundingBox &b = reference.bbox; + // Every extreme coordinate in either trajectory must be matched to a point + // in the other. This remains useful even when their boxes overlap. + const double dx = std::max(std::abs(a.x_min - b.x_min), std::abs(a.x_max - b.x_max)); + const double dy = std::max(std::abs(a.y_min - b.y_min), std::abs(a.y_max - b.y_max)); + if (!std::isfinite(dx) || !std::isfinite(dy)) + return 0.0; // Views without a cached box use unbounded defaults. + return std::max(dx * dx, dy * dy); +} + +double squared_points_to_box_lower_bound(const core::TrajectoryView &query, + const core::TrajectoryView &reference, double cutoff) { + if (query.size() == 0 || reference.size() == 0) { + return std::numeric_limits::infinity(); + } + double bound = 0.0; + for (std::size_t i = 0; i < query.size(); ++i) { + bound = std::max(bound, squared_dist_to_box(query.x[i], query.y[i], reference.bbox)); + if (bound >= cutoff) + return bound; + } + for (std::size_t j = 0; j < reference.size(); ++j) { + bound = std::max(bound, squared_dist_to_box(reference.x[j], reference.y[j], query.bbox)); + if (bound >= cutoff) + return bound; + } + return bound; +} + +// A greedy monotone path gives an upper bound, not an approximate answer. It +// supplies a finite pruning threshold before the exact top-k heap has filled. +double squared_greedy_upper_bound(const core::TrajectoryView &a, const core::TrajectoryView &b) { + if (a.size() == 0 || b.size() == 0) + return std::numeric_limits::infinity(); + const double inf = std::numeric_limits::infinity(); + std::size_t i = 0, j = 0; + double bound = squared_point_dist(a.x[0], a.y[0], b.x[0], b.y[0]); + while (i + 1 < a.size() || j + 1 < b.size()) { + const double diagonal = + i + 1 < a.size() && j + 1 < b.size() + ? squared_point_dist(a.x[i + 1], a.y[i + 1], b.x[j + 1], b.y[j + 1]) + : inf; + const double down = + i + 1 < a.size() ? squared_point_dist(a.x[i + 1], a.y[i + 1], b.x[j], b.y[j]) : inf; + const double right = + j + 1 < b.size() ? squared_point_dist(a.x[i], a.y[i], b.x[j + 1], b.y[j + 1]) : inf; + if (diagonal <= down && diagonal <= right) { + ++i; + ++j; + bound = std::max(bound, diagonal); + } else if (down <= right) { + ++i; + bound = std::max(bound, down); + } else { + ++j; + bound = std::max(bound, right); + } + } + return bound; +} + +struct Candidate { + double squared_lower_bound; + std::size_t idx; +}; + } // namespace struct neighbor { @@ -76,7 +181,8 @@ struct neighbor { double squared_frechet_distance; bool operator<(const neighbor &other) const { - return squared_frechet_distance < other.squared_frechet_distance; + return std::tie(squared_frechet_distance, idx) < + std::tie(other.squared_frechet_distance, other.idx); } }; @@ -86,30 +192,70 @@ 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. + // 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}; + } + 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); + }); + + double initial_cutoff = std::numeric_limits::infinity(); + if (k <= candidates.size()) { + constexpr std::size_t kInitialBoundCandidates = 64; + const auto count = std::min(candidates.size(), std::max(k, kInitialBoundCandidates)); + std::priority_queue upper_bounds; + for (std::size_t i = 0; i < count; ++i) { + const double upper = squared_greedy_upper_bound( + query_trajectory, trajectory_db.trajectory(candidates[i].idx)); + if (upper_bounds.size() < k) + upper_bounds.push(upper); + else if (upper < upper_bounds.top()) { + upper_bounds.pop(); + upper_bounds.push(upper); + } + } + // Include candidates equal to the upper bound, since we still need k results. + initial_cutoff = + std::nextafter(upper_bounds.top(), std::numeric_limits::infinity()); + } + std::priority_queue curr_neighbors; std::vector prev; std::vector curr; - for (std::size_t i = 0; i < trajectory_db.size(); ++i) { + for (const Candidate &candidate : candidates) { + const std::size_t i = candidate.idx; - core::TrajectoryView reference_trajectory = trajectory_db.trajectory(i); + // The greedy paths bound the k-th distance until exact results tighten it. + const double curr_max_dist = + curr_neighbors.size() < k + ? initial_cutoff + : std::nextafter(curr_neighbors.top().squared_frechet_distance, + std::numeric_limits::infinity()); + if (candidate.squared_lower_bound >= curr_max_dist) { + break; + } - // Until the heap is full any distance qualifies, so there is no bound to stop early at. - const double curr_max_dist = curr_neighbors.size() < k - ? std::numeric_limits::infinity() - : curr_neighbors.top().squared_frechet_distance; - // Even the best case for this trajectory cannot beat the current k-th neighbour. - // shortest leash can never be shorter than first or last point - if (curr_neighbors.size() == k && - squared_endpoint_lower_bound(query_trajectory, reference_trajectory) >= curr_max_dist) { + 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)) { continue; } double curr_dist = get_squared_frechet_distance(query_trajectory, reference_trajectory, curr_max_dist, prev, curr); + if (curr_dist >= curr_max_dist && std::isfinite(curr_max_dist)) + continue; if (curr_neighbors.size() < k) { curr_neighbors.push({.idx = i, .squared_frechet_distance = curr_dist}); continue; } - if (curr_dist < curr_max_dist) { + if (neighbor{.idx = i, .squared_frechet_distance = curr_dist} < curr_neighbors.top()) { curr_neighbors.pop(); curr_neighbors.push({.idx = i, .squared_frechet_distance = curr_dist}); } diff --git a/engine/tests/evaluator_test.cpp b/engine/tests/evaluator_test.cpp index 84080cb..8600451 100644 --- a/engine/tests/evaluator_test.cpp +++ b/engine/tests/evaluator_test.cpp @@ -100,8 +100,31 @@ TEST(Evaluator, CollectsKnnAndWhereAtScoresTogether) { const auto scores = evaluator.evaluate(); EXPECT_DOUBLE_EQ(scores.knn.kendalls_tau, 1); EXPECT_DOUBLE_EQ(scores.knn.recall_at_k, 1); + EXPECT_EQ(scores.knn.query_count, 2u); + EXPECT_EQ(scores.knn.tau_query_count, 2u); EXPECT_DOUBLE_EQ(scores.where_at.within_radius, 1); EXPECT_DOUBLE_EQ(scores.where_at.average_hit, 1); EXPECT_DOUBLE_EQ(evaluator.get_knn_scores().kendalls_tau, scores.knn.kendalls_tau); EXPECT_DOUBLE_EQ(evaluator.get_knn_scores().recall_at_k, scores.knn.recall_at_k); } + +TEST(KnnScores, ReportsCoverageWhenOnlySomeQueriesShareEnoughNeighbors) { + const core::TrajectoryDb db({0, 2, 4, 6, 8}, {10, 20, 30, 40}, {0, 10, 0, 10, 0, 10, 0, 10}, + std::vector(8, 39.95), + {116.35, 116.35, 116.36, 116.36, 116.37, 116.37, 116.38, 116.38}); + const query::KnnWorkload workload{db, {{0, {0, 1}}, {0, {0, 2}}, {0, {2, 3}}}, 2}; + const auto scores = evaluator::evaluate_knn(db, workload); + EXPECT_EQ(scores.query_count, 3u); + EXPECT_EQ(scores.tau_query_count, 1u); + EXPECT_DOUBLE_EQ(scores.kendalls_tau, 1.0); + EXPECT_DOUBLE_EQ(scores.recall_at_k, 0.5); +} + +TEST(KnnScores, EmptyWorkloadHasZeroCoverageAndUndefinedScores) { + const auto db = original(); + const auto scores = evaluator::evaluate_knn(db, {db, {}, 2}); + EXPECT_EQ(scores.query_count, 0u); + EXPECT_EQ(scores.tau_query_count, 0u); + EXPECT_TRUE(std::isnan(scores.kendalls_tau)); + EXPECT_TRUE(std::isnan(scores.recall_at_k)); +} diff --git a/engine/tests/knn_test.cpp b/engine/tests/knn_test.cpp index 6a68be9..4de26d5 100644 --- a/engine/tests/knn_test.cpp +++ b/engine/tests/knn_test.cpp @@ -33,6 +33,7 @@ core::TrajectoryDb make_db(const std::vector &tracks) { return {std::move(offsets), std::move(ids), std::move(t), std::move(lat), std::move(lon)}; } +// Random walks starting anywhere in a ~20 km area, so some are near each other and most are not std::vector random_tracks(std::mt19937_64 &rng, std::size_t count) { std::uniform_real_distribution start_lat(39.8, 40.0); std::uniform_real_distribution start_lon(116.2, 116.45); @@ -123,3 +124,39 @@ TEST(Knn, MatchesBruteForceOnRandomTrajectories) { } } } + +TEST(Knn, MatchesBruteForceForOverlappingLongTrajectories) { + std::mt19937_64 rng(19); + std::normal_distribution step(0.0, 0.0001); + std::vector tracks(28, Track(90)); + for (auto &track : tracks) { + double lat = 39.9, lon = 116.3; + for (auto &point : track) { + lat += step(rng); + lon += step(rng); + point = {lat, lon}; + } + } + const auto db = make_db({tracks.begin(), tracks.begin() + 20}); + const auto queries = make_db({tracks.begin() + 20, tracks.end()}); + for (const std::size_t k : {1u, 5u, 20u, 30u}) { + 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)); + auto query = db.trajectory(4); + query.bbox = {}; + EXPECT_EQ(query::knn(db, query, 5), brute_force_knn(db, query, 5)); +} + +TEST(Knn, BreaksEqualDistanceTiesByDatabaseIndex) { + const auto db = make_db(std::vector(80, {{39.9, 116.3}, {39.91, 116.31}})); + EXPECT_EQ(query::knn(db, db.trajectory(70), 5), (std::vector{0, 1, 2, 3, 4})); +} diff --git a/engine/tests/subset_test.cpp b/engine/tests/subset_test.cpp index 7254666..18fdc2d 100644 --- a/engine/tests/subset_test.cpp +++ b/engine/tests/subset_test.cpp @@ -39,6 +39,10 @@ void expect_same_trajectory(const core::TrajectoryView &a, const core::Trajector EXPECT_TRUE(std::ranges::equal(a.lon, b.lon)); EXPECT_TRUE(std::ranges::equal(a.x, b.x)); EXPECT_TRUE(std::ranges::equal(a.y, b.y)); + EXPECT_EQ(a.bbox.x_min, b.bbox.x_min); + EXPECT_EQ(a.bbox.x_max, b.bbox.x_max); + EXPECT_EQ(a.bbox.y_min, b.bbox.y_min); + EXPECT_EQ(a.bbox.y_max, b.bbox.y_max); } } // namespace @@ -98,3 +102,14 @@ TEST(Subset, RejectsOutOfRangeIndices) { const std::vector indices{1}; EXPECT_THROW(core::subset(db, indices), std::out_of_range); } + +TEST(TrajectoryDb, BoundingBoxCoversExactlyTheTrajectorysPoints) { + const auto db = make_db({{0, 1, 2}, {5, 6}}); + for (std::size_t i = 0; i < db.size(); ++i) { + const auto trajectory = db.trajectory(i); + EXPECT_EQ(trajectory.bbox.x_min, std::ranges::min(trajectory.x)); + EXPECT_EQ(trajectory.bbox.x_max, std::ranges::max(trajectory.x)); + EXPECT_EQ(trajectory.bbox.y_min, std::ranges::min(trajectory.y)); + EXPECT_EQ(trajectory.bbox.y_max, std::ranges::max(trajectory.y)); + } +}