From c8792e42b86a80a93a584c9bc792e315ee4ed2d6 Mon Sep 17 00:00:00 2001 From: Kristofferjoha Date: Sat, 7 Nov 2026 10:21:07 +0100 Subject: [PATCH 1/3] clustering --- README.md | 1 + engine/CMakeLists.txt | 3 + engine/configs/geolife.json | 7 +- .../backend/engine/evaluator/combined.hpp | 5 +- .../backend/engine/evaluator/evaluator.hpp | 5 + .../engine/evaluator/route_clustering.hpp | 34 +++ .../backend/engine/pipeline/config.hpp | 2 + .../backend/engine/query/route_clustering.hpp | 58 ++++ engine/src/evaluator/combined.cpp | 5 +- engine/src/evaluator/evaluator.cpp | 8 +- engine/src/evaluator/route_clustering.cpp | 95 +++++++ engine/src/pipeline/config.cpp | 16 +- engine/src/pipeline/pipeline.cpp | 30 +- engine/src/query/route_clustering.cpp | 124 +++++++++ engine/tests/config_test.cpp | 25 ++ engine/tests/route_clustering_test.cpp | 256 ++++++++++++++++++ 16 files changed, 665 insertions(+), 9 deletions(-) create mode 100644 engine/include/backend/engine/evaluator/route_clustering.hpp create mode 100644 engine/include/backend/engine/query/route_clustering.hpp create mode 100644 engine/src/evaluator/route_clustering.cpp create mode 100644 engine/src/query/route_clustering.cpp create mode 100644 engine/tests/route_clustering_test.cpp diff --git a/README.md b/README.md index c11673a..ea1058f 100644 --- a/README.md +++ b/README.md @@ -80,6 +80,7 @@ and it then falls back to its default. [engine/configs/geolife.json](engine/conf | `range` | Box size (m), time window (s) and sampling limits for range queries | | `where_at` | Probes per trajectory and max time gap (s) between points to probe | | `where_at_scoring` | Hit radius and tolerance range (m) for scoring where-at queries | +| `clustering` | Route clustering metric: `enabled`, neighbour threshold `eps_m` (m) and `min_size` per cluster. Reuses the self join when `join.mode` is `Self` with the same eps | | `allocation_alphas` | Allocation exponents to test: budgets follow size^α (`1` = size-proportional, `0` = equal). Exports for α ≠ 1 get an `_a<α>` suffix | ## Tests diff --git a/engine/CMakeLists.txt b/engine/CMakeLists.txt index a492dd9..5253712 100644 --- a/engine/CMakeLists.txt +++ b/engine/CMakeLists.txt @@ -37,6 +37,8 @@ add_library(engine src/query/join.cpp src/evaluator/join.cpp src/evaluator/combined.cpp + src/query/route_clustering.cpp + src/evaluator/route_clustering.cpp ) target_include_directories(engine PUBLIC include) @@ -69,6 +71,7 @@ add_executable(engine_tests tests/config_test.cpp tests/join_test.cpp tests/combined_test.cpp + tests/route_clustering_test.cpp ) target_link_libraries(engine_tests PRIVATE engine GTest::gtest_main Arrow::arrow_shared Parquet::parquet_shared) gtest_discover_tests(engine_tests) diff --git a/engine/configs/geolife.json b/engine/configs/geolife.json index 11fb1db..bf413f9 100644 --- a/engine/configs/geolife.json +++ b/engine/configs/geolife.json @@ -1,5 +1,5 @@ { - "data_path": "data/geolife.parquet", + "data_path": "data/geolife_gap60.parquet", "export_results": true, "export_dir": "data/output/geolife", "seed": 42, @@ -34,5 +34,10 @@ "join": { "mode": "Self", "eps_m": 500.0 + }, + "clustering": { + "enabled": false, + "eps_m": 150.0, + "min_size": 5 } } diff --git a/engine/include/backend/engine/evaluator/combined.hpp b/engine/include/backend/engine/evaluator/combined.hpp index 3d4256d..0472d46 100644 --- a/engine/include/backend/engine/evaluator/combined.hpp +++ b/engine/include/backend/engine/evaluator/combined.hpp @@ -7,7 +7,8 @@ namespace evaluator { // 0 = the floor (endpoints only), 1 = perfect, without clipping. Combined averages -// one primary score per type: recall@k, range F1, where-at within-radius hit rate, join F1. +// one primary score per type: recall@k, range F1, where-at within-radius hit rate, join F1 +// and route clustering. // Missing primaries are skipped, never replaced with a diagnostic score. // The radius is configured by hit_radius_m (100 m for the research protocol). // Average-hit remains a diagnostic and does not contribute to the combined score. @@ -21,11 +22,13 @@ struct NormalizedScores { double where_at_within_radius = missing; double where_at_average_hit = missing; double join_f1 = missing; + double clustering_score = missing; double knn = missing; double range = missing; double where_at = missing; double join = missing; + double clustering = missing; double combined = missing; }; diff --git a/engine/include/backend/engine/evaluator/evaluator.hpp b/engine/include/backend/engine/evaluator/evaluator.hpp index dbb2d70..f2f38f0 100644 --- a/engine/include/backend/engine/evaluator/evaluator.hpp +++ b/engine/include/backend/engine/evaluator/evaluator.hpp @@ -3,10 +3,12 @@ #include "backend/engine/core/trajectory_db.hpp" #include "backend/engine/evaluator/join.hpp" #include "backend/engine/evaluator/knn.hpp" +#include "backend/engine/evaluator/route_clustering.hpp" #include "backend/engine/evaluator/where_at.hpp" #include "backend/engine/query/join.hpp" #include "backend/engine/query/knn.hpp" #include "backend/engine/query/range_query.hpp" +#include "backend/engine/query/route_clustering.hpp" #include #include @@ -19,6 +21,7 @@ struct Workloads { std::span where_at = {}; std::span range = {}; const query::JoinWorkload *join = nullptr; // a pointer, as it also carries eps and the queries + const query::RouteClusteringWorkload *clustering = nullptr; }; // Kendall's tau is in [-1, 1] and the others in [0, 1], or NaN if it could not be computed @@ -28,6 +31,7 @@ struct Scoretable { WhereAtScores where_at; JoinScores join; double range_f1 = std::numeric_limits::quiet_NaN(); + RouteClusteringScores clustering; }; struct EvaluationResult { @@ -49,6 +53,7 @@ class Evaluator { double get_range_f1() const; WhereAtScores get_where_at_scores() const; JoinScores get_join_scores() const; + RouteClusteringScores get_clustering_scores() const; private: const core::TrajectoryDb &m_simplified; diff --git a/engine/include/backend/engine/evaluator/route_clustering.hpp b/engine/include/backend/engine/evaluator/route_clustering.hpp new file mode 100644 index 0000000..940f238 --- /dev/null +++ b/engine/include/backend/engine/evaluator/route_clustering.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include "backend/engine/core/trajectory_db.hpp" +#include "backend/engine/query/route_clustering.hpp" + +#include +#include + +namespace evaluator { + +// Reruns the route clustering on the compressed database and compares it with the reference. +// score is the symmetric macro best-match F1: each cluster's best overlap F1 with any cluster on +// the other side, averaged per side, then over both sides. NaN without reference clusters, 0 if +// nothing is predicted. +struct RouteClusteringScores { + static constexpr double missing = std::numeric_limits::quiet_NaN(); + + double score = missing; + std::size_t reference_clusters = 0; + std::size_t predicted_clusters = 0; + std::size_t reference_covered = 0; + std::size_t predicted_covered = 0; +}; + +query::RouteClusters compressed_route_clusters(const core::TrajectoryDb &compressed, + const query::RouteClusteringWorkload &workload); + +RouteClusteringScores compare_route_clusters(const query::RouteClusters &reference, + const query::RouteClusters &predicted); + +RouteClusteringScores evaluate_route_clustering(const core::TrajectoryDb &compressed, + const query::RouteClusteringWorkload &workload); + +} // namespace evaluator diff --git a/engine/include/backend/engine/pipeline/config.hpp b/engine/include/backend/engine/pipeline/config.hpp index 965903f..74f47db 100644 --- a/engine/include/backend/engine/pipeline/config.hpp +++ b/engine/include/backend/engine/pipeline/config.hpp @@ -4,6 +4,7 @@ #include "backend/engine/query/join.hpp" #include "backend/engine/query/knn.hpp" #include "backend/engine/query/range_query.hpp" +#include "backend/engine/query/route_clustering.hpp" #include "backend/engine/query/where_at_query.hpp" #include "backend/engine/selection/selector.hpp" #include "backend/engine/split/split.hpp" @@ -29,6 +30,7 @@ struct Config { query::WhereAtConfig where_at; evaluator::WhereAtScoringConfig where_at_scoring; query::JoinConfig join; + query::RouteClusteringConfig clustering; }; // Reads a JSON config. Only data_path is required, a missing key keeps its default. diff --git a/engine/include/backend/engine/query/route_clustering.hpp b/engine/include/backend/engine/query/route_clustering.hpp new file mode 100644 index 0000000..2161372 --- /dev/null +++ b/engine/include/backend/engine/query/route_clustering.hpp @@ -0,0 +1,58 @@ +#pragma once + +#include "backend/engine/core/trajectory_db.hpp" +#include "backend/engine/query/join.hpp" + +#include +#include +#include +#include +#include + +namespace query { + +struct RouteClusteringConfig { + bool enabled = false; + double eps_m = 500.0; // neighbour threshold, same test as the self join + std::size_t min_size = 5; // a cluster needs at least this many trajectories, medoid included +}; + +using TrajectoryPair = std::pair; + +// Greedy medoid grouping. Repeatedly takes the unassigned trajectory with the most unassigned +// neighbours and groups it with them if that gives at least min_size trajectories. Every member +// neighbours its medoid, so clusters cannot chain. Ties go to the earliest start time, then the +// smallest start x, then the smallest trajectory ID, which makes the result order independent. +struct RouteClusters { + std::vector label; // cluster per trajectory, -1 = unclustered + std::vector> members; // members[c][0] is the medoid +}; + +RouteClusters cluster_routes(const core::TrajectoryDb &trajectories, + std::span neighbours, std::size_t min_size); + +// ID and endpoints of an original trajectory. The candidate pairs were pruned by endpoint +// distance, so a database is only comparable if every one of these matches. +struct TrajectoryEndpoints { + std::int32_t trajectory_id = 0; + double first_x = 0.0, first_y = 0.0, last_x = 0.0, last_y = 0.0; + bool operator==(const TrajectoryEndpoints &) const = default; +}; + +TrajectoryEndpoints endpoints_of(const core::TrajectoryView &trajectory); + +struct RouteClusteringWorkload { + std::vector candidates; // endpoint-pruned pairs, i < j + RouteClusters reference; + std::vector originals; + double eps_m = 0.0; + std::size_t min_size = 0; +}; + +// Clusters the original database. A self join workload over the same trajectories and eps is +// reused instead of recomputing the Fréchet tests on the originals. +RouteClusteringWorkload generate_route_clustering(const core::TrajectoryDb &original, + const RouteClusteringConfig &config = {}, + const JoinWorkload *self_join = nullptr); + +} // namespace query diff --git a/engine/src/evaluator/combined.cpp b/engine/src/evaluator/combined.cpp index 2be9b7b..19a8f27 100644 --- a/engine/src/evaluator/combined.cpp +++ b/engine/src/evaluator/combined.cpp @@ -41,13 +41,16 @@ NormalizedScores normalize_scores(const Scoretable &scores, const Scoretable &fl result.where_at_average_hit = normalize(scores.where_at.average_hit, floor.where_at.average_hit); result.join_f1 = normalize(scores.join.f1, floor.join.f1); + result.clustering_score = normalize(scores.clustering.score, floor.clustering.score); result.knn = result.knn_recall_at_k; result.range = result.range_f1; result.where_at = result.where_at_within_radius; result.join = result.join_f1; + result.clustering = result.clustering_score; - result.combined = mean_ignoring_nan({result.knn, result.range, result.where_at, result.join}); + result.combined = mean_ignoring_nan( + {result.knn, result.range, result.where_at, result.join, result.clustering}); return result; } diff --git a/engine/src/evaluator/evaluator.cpp b/engine/src/evaluator/evaluator.cpp index 9faf8a7..5528124 100644 --- a/engine/src/evaluator/evaluator.cpp +++ b/engine/src/evaluator/evaluator.cpp @@ -27,7 +27,8 @@ EvaluationResult Evaluator::evaluate() const { .range_recall = range.recall, .where_at = get_where_at_scores(), .join = get_join_scores(), - .range_f1 = range.f1}, + .range_f1 = range.f1, + .clustering = get_clustering_scores()}, .knn = std::move(knn.answers), .range = std::move(range.answers)}; } @@ -52,4 +53,9 @@ JoinScores Evaluator::get_join_scores() const { return m_workloads.join ? evaluate_join(m_simplified, *m_workloads.join) : JoinScores{}; } +RouteClusteringScores Evaluator::get_clustering_scores() const { + return m_workloads.clustering ? evaluate_route_clustering(m_simplified, *m_workloads.clustering) + : RouteClusteringScores{}; +} + } // namespace evaluator diff --git a/engine/src/evaluator/route_clustering.cpp b/engine/src/evaluator/route_clustering.cpp new file mode 100644 index 0000000..4a30cd9 --- /dev/null +++ b/engine/src/evaluator/route_clustering.cpp @@ -0,0 +1,95 @@ +#include "backend/engine/evaluator/route_clustering.hpp" +#include "backend/engine/core/parallel.hpp" +#include "backend/engine/query/frechet_helper.hpp" + +#include +#include +#include +#include +#include + +namespace evaluator { + +namespace { + +// Mean over a's clusters of the best overlap F1 with any of b's clusters +double best_match(const query::RouteClusters &a, const query::RouteClusters &b) { + double sum = 0.0; + for (const auto &members : a.members) { + std::unordered_map overlap; + for (const auto i : members) { + if (b.label[i] >= 0) + ++overlap[b.label[i]]; + } + double best = 0.0; + for (const auto &[cluster, shared] : overlap) { + const auto size = members.size() + b.members[cluster].size(); + best = std::max(best, 2.0 * static_cast(shared) / static_cast(size)); + } + sum += best; + } + return sum / static_cast(a.members.size()); +} + +std::size_t covered(const query::RouteClusters &clusters) { + std::size_t total = 0; + for (const auto &members : clusters.members) + total += members.size(); + return total; +} + +} // namespace + +query::RouteClusters compressed_route_clusters(const core::TrajectoryDb &compressed, + const query::RouteClusteringWorkload &workload) { + if (compressed.size() != workload.originals.size()) + throw std::invalid_argument("route clustering: database size does not match reference"); + // Simplification keeps endpoints, so the original candidates cover every possible match + for (std::size_t i = 0; i < compressed.size(); ++i) { + const auto trajectory = compressed.trajectory(i); + if (trajectory.size() == 0 || query::endpoints_of(trajectory) != workload.originals[i]) + throw std::invalid_argument( + "route clustering: IDs/order or endpoints do not match reference"); + } + + std::vector matched(workload.candidates.size(), 0); + core::parallel_for(workload.candidates.size(), [&](std::size_t k) { + const auto [i, j] = workload.candidates[k]; + matched[k] = query::frechet_within(compressed.trajectory(i), compressed.trajectory(j), + workload.eps_m); + }); + std::vector neighbours; + for (std::size_t k = 0; k < matched.size(); ++k) { + if (matched[k]) + neighbours.push_back(workload.candidates[k]); + } + return query::cluster_routes(compressed, neighbours, workload.min_size); +} + +RouteClusteringScores compare_route_clusters(const query::RouteClusters &reference, + const query::RouteClusters &predicted) { + if (reference.label.size() != predicted.label.size()) + throw std::invalid_argument("route clustering: label sizes differ"); + + RouteClusteringScores scores; + scores.reference_clusters = reference.members.size(); + scores.predicted_clusters = predicted.members.size(); + scores.reference_covered = covered(reference); + scores.predicted_covered = covered(predicted); + if (reference.members.empty()) + return scores; + + scores.score = + predicted.members.empty() + ? 0.0 + : 0.5 * (best_match(reference, predicted) + best_match(predicted, reference)); + return scores; +} + +RouteClusteringScores evaluate_route_clustering(const core::TrajectoryDb &compressed, + const query::RouteClusteringWorkload &workload) { + return compare_route_clusters(workload.reference, + compressed_route_clusters(compressed, workload)); +} + +} // namespace evaluator diff --git a/engine/src/pipeline/config.cpp b/engine/src/pipeline/config.cpp index bc48c14..0785f64 100644 --- a/engine/src/pipeline/config.cpp +++ b/engine/src/pipeline/config.cpp @@ -92,11 +92,24 @@ void parse_join_config(Config &config, const json &join) { config.join.eps_m = join.value("eps_m", config.join.eps_m); } +void parse_clustering_config(Config &config, const json &clustering) { + validate_keys(clustering, {"enabled", "eps_m", "min_size"}, "clustering"); + config.clustering.enabled = clustering.value("enabled", config.clustering.enabled); + config.clustering.eps_m = clustering.value("eps_m", config.clustering.eps_m); + if (clustering.contains("min_size")) { + // Checked before conversion, since a negative value would wrap around in size_t + const auto &min_size = clustering["min_size"]; + if (!min_size.is_number_unsigned() || min_size.get() < 2) + throw std::invalid_argument("config: clustering.min_size must be an integer >= 2"); + config.clustering.min_size = min_size.get(); + } +} + Config parse_config(const json &data) { validate_keys(data, {"data_path", "export_results", "export_dir", "seed", "budget_shares", "allocation_alphas", "algorithms", "split", "knn", "range", "where_at", - "where_at_scoring", "join"}, + "where_at_scoring", "join", "clustering"}, "root"); Config config; @@ -130,6 +143,7 @@ Config parse_config(const json &data) { parse_where_at_config(config, data.value("where_at", json::object())); parse_where_at_scoring_config(config, data.value("where_at_scoring", json::object())); parse_join_config(config, data.value("join", json::object())); + parse_clustering_config(config, data.value("clustering", json::object())); return config; } diff --git a/engine/src/pipeline/pipeline.cpp b/engine/src/pipeline/pipeline.cpp index b566f9d..3ffe51c 100644 --- a/engine/src/pipeline/pipeline.cpp +++ b/engine/src/pipeline/pipeline.cpp @@ -9,6 +9,7 @@ #include "backend/engine/query/join.hpp" #include "backend/engine/query/knn.hpp" #include "backend/engine/query/range_query.hpp" +#include "backend/engine/query/route_clustering.hpp" #include "backend/engine/query/where_at_query.hpp" #include "backend/engine/selection/selector.hpp" #include "backend/engine/split/split.hpp" @@ -59,12 +60,18 @@ void print_scores(const evaluator::Scoretable &scores) { std::cout << " range F1: " << scores.range_f1 << ", recall: " << scores.range_recall << '\n'; std::cout << " join F1: " << scores.join.f1 << ", precision: " << scores.join.precision << ", recall: " << scores.join.recall << '\n'; + if (!std::isnan(scores.clustering.score) || scores.clustering.reference_clusters != 0) + std::cout << " route clustering: " << scores.clustering.score << " (" + << scores.clustering.predicted_clusters << " clusters, " + << scores.clustering.predicted_covered << " trajectories)\n"; } void print_normalized(const evaluator::NormalizedScores &scores) { std::cout << " normalized: knn=" << scores.knn << ", range=" << scores.range - << ", where-at=" << scores.where_at << ", join=" << scores.join - << ", combined=" << scores.combined << '\n'; + << ", where-at=" << scores.where_at << ", join=" << scores.join; + if (!std::isnan(scores.clustering)) + std::cout << ", clustering=" << scores.clustering; + std::cout << ", combined=" << scores.combined << '\n'; } core::Simplification endpoints_only(const core::TrajectoryDb &trajectories) { @@ -105,13 +112,28 @@ void run(const Config &config) { : " queries, eps=") << join.eps_m << " m, " << join_pairs << " matching pairs\n"; + query::RouteClusteringWorkload clustering; + if (config.clustering.enabled) { + clustering = query::generate_route_clustering(database, config.clustering, &join); + std::size_t covered = 0; + for (const auto &members : clustering.reference.members) + covered += members.size(); + std::cout << "route clustering: " << clustering.reference.members.size() + << " clusters covering " << covered << " trajectories, eps=" << clustering.eps_m + << " m, min size " << clustering.min_size << '\n'; + } + if (config.export_results) { io::export_workloads(config.export_dir, database, range.queries, knn_workload); std::cout << "workloads written to " << config.export_dir.string() << '\n'; } - const evaluator::Workloads workloads{ - .knn = &knn_workload, .where_at = where_at.queries, .range = range.queries, .join = &join}; + const evaluator::Workloads workloads{.knn = &knn_workload, + .where_at = where_at.queries, + .range = range.queries, + .join = &join, + .clustering = + config.clustering.enabled ? &clustering : nullptr}; const auto floor_db = core::materialize(database, endpoints_only(database)); std::cout << "floor, endpoints only: "; diff --git a/engine/src/query/route_clustering.cpp b/engine/src/query/route_clustering.cpp new file mode 100644 index 0000000..b173b8f --- /dev/null +++ b/engine/src/query/route_clustering.cpp @@ -0,0 +1,124 @@ +#include "backend/engine/query/route_clustering.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace query { + +TrajectoryEndpoints endpoints_of(const core::TrajectoryView &trajectory) { + if (trajectory.size() == 0) + throw std::invalid_argument("route clustering: empty trajectory"); + return {trajectory.trajectory_id, trajectory.x.front(), trajectory.y.front(), + trajectory.x.back(), trajectory.y.back()}; +} + +RouteClusters cluster_routes(const core::TrajectoryDb &trajectories, + std::span neighbours, std::size_t min_size) { + if (min_size < 2) + throw std::invalid_argument("route clustering: min_size must be at least 2"); + const std::size_t n = trajectories.size(); + std::vector> adjacency(n); + for (const auto &[i, j] : neighbours) { + if (i >= n || j >= n || i == j) + throw std::invalid_argument("route clustering: invalid neighbour pair"); + adjacency[i].push_back(j); + adjacency[j].push_back(i); + } + // A repeated pair must not count twice towards the degree + for (auto &list : adjacency) { + std::ranges::sort(list); + list.erase(std::ranges::unique(list).begin(), list.end()); + } + + // Max-heap on (degree, then tie-breaks). Degrees only fall, so stale entries are re-pushed. + using Key = std::tuple; + const auto key = [&](std::size_t i, std::size_t degree) { + const auto trajectory = trajectories.trajectory(i); + if (trajectory.size() == 0) + throw std::invalid_argument("route clustering: empty trajectory"); + return Key{degree, -trajectory.t.front(), -trajectory.x.front(), + -static_cast(trajectory.trajectory_id), i}; + }; + std::priority_queue heap; + for (std::size_t i = 0; i < n; ++i) + heap.push(key(i, adjacency[i].size())); + + RouteClusters clusters; + clusters.label.assign(n, -1); + while (!heap.empty()) { + const auto [degree, neg_t, neg_x, neg_id, i] = heap.top(); + heap.pop(); + if (clusters.label[i] >= 0) + continue; + std::size_t free = 0; + for (const auto j : adjacency[i]) + free += clusters.label[j] < 0; + if (free < degree) { + heap.push(key(i, free)); + continue; + } + if (free + 1 < min_size) + break; + const auto id = static_cast(clusters.members.size()); + auto &members = clusters.members.emplace_back(); + members.push_back(i); + clusters.label[i] = id; + for (const auto j : adjacency[i]) { + if (clusters.label[j] < 0) { + clusters.label[j] = id; + members.push_back(j); + } + } + } + return clusters; +} + +RouteClusteringWorkload generate_route_clustering(const core::TrajectoryDb &original, + const RouteClusteringConfig &config, + const JoinWorkload *self_join) { + if (!std::isfinite(config.eps_m) || config.eps_m <= 0.0) + throw std::invalid_argument("route clustering: eps_m must be finite and positive"); + + RouteClusteringWorkload workload; + workload.eps_m = config.eps_m; + workload.min_size = config.min_size; + workload.originals.reserve(original.size()); + for (std::size_t i = 0; i < original.size(); ++i) + workload.originals.push_back(endpoints_of(original.trajectory(i))); + + // A self join stores its database's IDs and endpoints, which must be these trajectories' + const auto same_database = [&] { + if (self_join->queries.size() != original.size() || + self_join->entries.size() != original.size()) + return false; + for (std::size_t i = 0; i < original.size(); ++i) { + const auto stored = self_join->queries.trajectory(i); + if (stored.size() == 0 || endpoints_of(stored) != workload.originals[i]) + return false; + } + return true; + }; + std::vector reference; + const bool reuse = self_join && self_join->mode == JoinMode::Self && + self_join->eps_m == config.eps_m && same_database(); + if (reuse) { + for (const auto &entry : self_join->entries) { + for (const auto j : entry.candidates) + workload.candidates.emplace_back(entry.query_index, j); + for (const auto j : entry.reference) + reference.emplace_back(entry.query_index, j); + } + } else { + const auto join = generate_self_join(original, config.eps_m); + return generate_route_clustering(original, config, &join); + } + workload.reference = cluster_routes(original, reference, config.min_size); + return workload; +} + +} // namespace query diff --git a/engine/tests/config_test.cpp b/engine/tests/config_test.cpp index b6cc8cb..c8f31c3 100644 --- a/engine/tests/config_test.cpp +++ b/engine/tests/config_test.cpp @@ -182,3 +182,28 @@ TEST(Config, RejectsEmptyAllocationAlphas) { "no_alphas.json", R"({"data_path": "d", "allocation_alphas": []})")), std::invalid_argument); } + +TEST(Config, ReadsClusteringKeys) { + const auto config = pipeline::load_config(write_config("clustering.json", R"({ + "data_path": "d", "clustering": {"enabled": true, "eps_m": 250.0, "min_size": 4} + })")); + EXPECT_TRUE(config.clustering.enabled); + EXPECT_DOUBLE_EQ(config.clustering.eps_m, 250.0); + EXPECT_EQ(config.clustering.min_size, 4u); +} + +TEST(Config, RejectsInvalidClustering) { + EXPECT_THROW( + pipeline::load_config(write_config("clustering_small.json", + R"({"data_path": "d", "clustering": {"min_size": 1}})")), + std::invalid_argument); + EXPECT_THROW(pipeline::load_config(write_config( + "clustering_key.json", R"({"data_path": "d", "clustering": {"minsize": 5}})")), + std::invalid_argument); + for (const char *value : {"-1", "2.5", "\"5\""}) + EXPECT_THROW( + pipeline::load_config(write_config( + "clustering_bad.json", + std::string(R"({"data_path": "d", "clustering": {"min_size": )") + value + "}}")), + std::invalid_argument); +} diff --git a/engine/tests/route_clustering_test.cpp b/engine/tests/route_clustering_test.cpp new file mode 100644 index 0000000..2616e86 --- /dev/null +++ b/engine/tests/route_clustering_test.cpp @@ -0,0 +1,256 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "backend/engine/core/simplification.hpp" +#include "backend/engine/core/trajectory_db.hpp" +#include "backend/engine/evaluator/combined.hpp" +#include "backend/engine/evaluator/route_clustering.hpp" +#include "backend/engine/query/join.hpp" +#include "backend/engine/query/route_clustering.hpp" + +namespace { + +using Track = std::vector>; // (lat, lon) per point + +core::TrajectoryDb make_db(const std::vector &tracks, + const std::vector &ids = {}) { + std::vector offsets{0}; + std::vector trajectory_ids; + std::vector t; + std::vector lat; + std::vector lon; + for (std::size_t i = 0; i < tracks.size(); ++i) { + trajectory_ids.push_back(ids.empty() ? static_cast(i) : ids[i]); + for (std::size_t k = 0; k < tracks[i].size(); ++k) { + t.push_back(static_cast(k)); + lat.push_back(tracks[i][k].first); + lon.push_back(tracks[i][k].second); + } + offsets.push_back(t.size()); + } + return {std::move(offsets), std::move(trajectory_ids), std::move(t), std::move(lat), + std::move(lon)}; +} + +// One short track per trajectory, start points spread along longitude so tie-breaks are distinct +core::TrajectoryDb spread_db(std::size_t n) { + std::vector tracks; + for (std::size_t i = 0; i < n; ++i) { + const double lon = 116.30 + 0.01 * static_cast(i); + tracks.push_back({{39.90, lon}, {39.91, lon}}); + } + return make_db(tracks); +} + +std::set> id_sets(const core::TrajectoryDb &db, + const query::RouteClusters &clusters) { + std::set> result; + for (const auto &members : clusters.members) { + std::set ids; + for (const auto i : members) + ids.insert(db.trajectory(i).trajectory_id); + result.insert(ids); + } + return result; +} + +// Trajectories on a few shared routes with small offsets, so routes form clusters +std::vector route_tracks(std::mt19937_64 &rng, std::size_t routes, std::size_t per_route) { + std::normal_distribution jitter(0.0, 0.0003); + std::vector tracks; + for (std::size_t r = 0; r < routes; ++r) { + for (std::size_t m = 0; m < per_route; ++m) { + Track track; + for (std::size_t k = 0; k < 20; ++k) { + const double s = static_cast(k) / 19.0; + track.emplace_back(39.90 + 0.05 * static_cast(r) + 0.02 * s + jitter(rng), + 116.30 + 0.03 * s + 0.01 * std::sin(6.0 * s) + jitter(rng)); + } + tracks.push_back(track); + } + } + return tracks; +} + +} // namespace + +TEST(RouteClustering, StarFormsOneClusterAndLeavesNoiseUnclustered) { + const auto db = spread_db(7); + const std::vector pairs{{0, 1}, {0, 2}, {0, 3}, {0, 4}}; + const auto clusters = query::cluster_routes(db, pairs, 5); + ASSERT_EQ(clusters.members.size(), 1u); + EXPECT_EQ(clusters.members[0].front(), 0u); + EXPECT_EQ(std::set(clusters.members[0].begin(), clusters.members[0].end()), + (std::set{0, 1, 2, 3, 4})); + EXPECT_EQ(clusters.label[5], -1); + EXPECT_EQ(clusters.label[6], -1); + EXPECT_TRUE(query::cluster_routes(db, pairs, 6).members.empty()); +} + +TEST(RouteClustering, MembersMustNeighbourTheMedoidSoPathsDoNotChain) { + const auto db = spread_db(5); + const std::vector path{{0, 1}, {1, 2}, {2, 3}, {3, 4}}; + const auto clusters = query::cluster_routes(db, path, 3); + for (const auto &members : clusters.members) { + EXPECT_LE(members.size(), 3u); + const auto medoid = members.front(); + for (std::size_t k = 1; k < members.size(); ++k) { + const auto a = std::min(medoid, members[k]); + const auto b = std::max(medoid, members[k]); + EXPECT_NE(std::find(path.begin(), path.end(), query::TrajectoryPair{a, b}), path.end()); + } + } + // Ties on degree go to the smallest start x, which is trajectory 1 + ASSERT_EQ(clusters.members.size(), 1u); + EXPECT_EQ(clusters.members[0].front(), 1u); +} + +TEST(RouteClustering, ResultDoesNotDependOnTrajectoryOrder) { + std::mt19937_64 rng(7); + const auto tracks = route_tracks(rng, 4, 6); + std::vector ids(tracks.size()); + std::iota(ids.begin(), ids.end(), 0); + const auto db = make_db(tracks, ids); + + std::vector order(tracks.size()); + std::iota(order.begin(), order.end(), 0); + std::ranges::shuffle(order, rng); + std::vector shuffled; + std::vector shuffled_ids; + for (const auto i : order) { + shuffled.push_back(tracks[i]); + shuffled_ids.push_back(ids[i]); + } + const auto permuted = make_db(shuffled, shuffled_ids); + + const auto a = query::generate_route_clustering(db, {.enabled = true, .min_size = 3}); + const auto b = query::generate_route_clustering(permuted, {.enabled = true, .min_size = 3}); + EXPECT_FALSE(a.reference.members.empty()); + EXPECT_EQ(id_sets(db, a.reference), id_sets(permuted, b.reference)); +} + +TEST(RouteClustering, ReusesAMatchingSelfJoin) { + std::mt19937_64 rng(3); + const auto db = make_db(route_tracks(rng, 3, 5)); + const auto join = query::generate_self_join(db, 500.0); + const auto reused = query::generate_route_clustering( + db, {.enabled = true, .eps_m = 500.0, .min_size = 3}, &join); + const auto fresh = + query::generate_route_clustering(db, {.enabled = true, .eps_m = 500.0, .min_size = 3}); + EXPECT_EQ(reused.candidates, fresh.candidates); + EXPECT_EQ(reused.reference.members, fresh.reference.members); +} + +TEST(RouteClustering, IgnoresASelfJoinOfAnotherDatabaseWithTheSameSize) { + std::mt19937_64 rng(3); + const auto db = make_db(route_tracks(rng, 3, 5)); + const auto other = make_db(std::vector(db.size(), {{39.0, 116.0}, {39.0, 116.001}})); + const auto other_join = query::generate_self_join(other, 500.0); + const query::RouteClusteringConfig config{.enabled = true, .eps_m = 500.0, .min_size = 3}; + const auto workload = query::generate_route_clustering(db, config, &other_join); + const auto fresh = query::generate_route_clustering(db, config); + EXPECT_EQ(workload.candidates, fresh.candidates); + EXPECT_EQ(workload.reference.members, fresh.reference.members); +} + +TEST(RouteClustering, RepeatedPairsCountOnce) { + const auto db = spread_db(3); + const std::vector repeated{{0, 1}, {1, 0}, {0, 1}}; + EXPECT_TRUE(query::cluster_routes(db, repeated, 3).members.empty()); + const std::vector star{{0, 1}, {1, 0}, {0, 2}}; + const auto clusters = query::cluster_routes(db, star, 3); + ASSERT_EQ(clusters.members.size(), 1u); + EXPECT_EQ(clusters.members[0].size(), 3u); +} + +TEST(RouteClustering, UncompressedDatabaseScoresPerfectly) { + std::mt19937_64 rng(11); + const auto db = make_db(route_tracks(rng, 4, 6)); + const auto workload = query::generate_route_clustering(db, {.enabled = true, .min_size = 3}); + ASSERT_FALSE(workload.reference.members.empty()); + const auto scores = evaluator::evaluate_route_clustering(db, workload); + EXPECT_DOUBLE_EQ(scores.score, 1.0); + EXPECT_EQ(scores.predicted_clusters, scores.reference_clusters); +} + +TEST(RouteClustering, CompressedDatabaseIsScoredAgainstTheReference) { + std::mt19937_64 rng(5); + const auto db = make_db(route_tracks(rng, 4, 6)); + const auto workload = query::generate_route_clustering(db, {.enabled = true, .min_size = 3}); + core::Simplification endpoints; + for (std::size_t i = 0; i < db.size(); ++i) { + const std::vector kept{0, db.trajectory(i).size() - 1}; + endpoints.add_trajectory(kept); + } + const auto compressed = core::materialize(db, endpoints); + const auto scores = evaluator::evaluate_route_clustering(compressed, workload); + EXPECT_GE(scores.score, 0.0); + EXPECT_LE(scores.score, 1.0); + const auto direct = evaluator::compare_route_clusters( + workload.reference, evaluator::compressed_route_clusters(compressed, workload)); + EXPECT_DOUBLE_EQ(scores.score, direct.score); +} + +TEST(RouteClustering, CompareScoresOverlapFromBothSides) { + query::RouteClusters reference{{0, 0, 0, 0, -1}, {{0, 1, 2, 3}}}; + query::RouteClusters predicted{{0, 0, -1, -1, -1}, {{0, 1}}}; + const auto scores = evaluator::compare_route_clusters(reference, predicted); + // Both sides' best F1 is 2 * 2 / (4 + 2) + EXPECT_NEAR(scores.score, 2.0 / 3.0, 1e-12); + EXPECT_DOUBLE_EQ(evaluator::compare_route_clusters(reference, reference).score, 1.0); + + const query::RouteClusters none{{-1, -1, -1, -1, -1}, {}}; + EXPECT_DOUBLE_EQ(evaluator::compare_route_clusters(reference, none).score, 0.0); + EXPECT_TRUE(std::isnan(evaluator::compare_route_clusters(none, predicted).score)); +} + +TEST(RouteClustering, RejectsMismatchedDatabasesAndInvalidInput) { + std::mt19937_64 rng(2); + const auto db = make_db(route_tracks(rng, 2, 4)); + const auto workload = query::generate_route_clustering(db, {.enabled = true, .min_size = 3}); + const auto other = make_db(route_tracks(rng, 2, 3)); + EXPECT_THROW(evaluator::evaluate_route_clustering(other, workload), std::invalid_argument); + // Same IDs and order, but a moved endpoint could match outside the cached candidates + std::vector endpoints; + for (std::size_t i = 0; i < db.size(); ++i) { + const auto t = db.trajectory(i); + endpoints.push_back({{t.lat.front(), t.lon.front()}, {t.lat.back(), t.lon.back()}}); + } + EXPECT_NO_THROW(evaluator::evaluate_route_clustering(make_db(endpoints), workload)); + endpoints[0].back().second += 0.01; + EXPECT_THROW(evaluator::evaluate_route_clustering(make_db(endpoints), workload), + std::invalid_argument); + EXPECT_THROW(query::cluster_routes(db, {}, 1), std::invalid_argument); + const std::vector self_pair{{0, 0}}; + EXPECT_THROW(query::cluster_routes(db, self_pair, 3), std::invalid_argument); + EXPECT_THROW(query::generate_route_clustering(db, {.enabled = true, .eps_m = -1.0}), + std::invalid_argument); +} + +TEST(RouteClustering, CombinedIncludesClustering) { + evaluator::Scoretable floor{}; + floor.knn.recall_at_k = 0.0; + floor.range_f1 = 0.0; + floor.where_at.within_radius = 0.0; + floor.join.f1 = 0.0; + floor.clustering.score = 0.5; + auto scores = floor; + scores.knn.recall_at_k = 1.0; + scores.clustering.score = 0.75; + const auto normalized = evaluator::normalize_scores(scores, floor); + EXPECT_DOUBLE_EQ(normalized.clustering, 0.5); + EXPECT_DOUBLE_EQ(normalized.combined, 1.5 / 5.0); + scores.clustering.score = std::numeric_limits::quiet_NaN(); + EXPECT_DOUBLE_EQ(evaluator::normalize_scores(scores, floor).combined, 0.25); +} From 4b1c79962c9e6119ae808d71ec9ef3b16db67375 Mon Sep 17 00:00:00 2001 From: Kristofferjoha Date: Sat, 7 Nov 2026 10:50:41 +0100 Subject: [PATCH 2/3] frechet distance --- README.md | 1 + engine/CMakeLists.txt | 2 + engine/configs/geolife.json | 7 +- .../backend/engine/query/frechet_helper.hpp | 15 + engine/include/backend/engine/query/join.hpp | 9 +- engine/include/backend/engine/query/knn.hpp | 11 +- .../backend/engine/query/route_clustering.hpp | 6 +- engine/src/evaluator/join.cpp | 3 +- engine/src/evaluator/knn.cpp | 2 +- engine/src/evaluator/route_clustering.cpp | 2 +- engine/src/pipeline/config.cpp | 9 +- engine/src/query/continuous_frechet.cpp | 200 ++++++++++++ engine/src/query/join.cpp | 19 +- engine/src/query/knn.cpp | 299 +++++++++++++---- engine/src/query/route_clustering.cpp | 6 +- engine/tests/config_test.cpp | 14 + engine/tests/continuous_frechet_test.cpp | 307 ++++++++++++++++++ 17 files changed, 828 insertions(+), 84 deletions(-) create mode 100644 engine/src/query/continuous_frechet.cpp create mode 100644 engine/tests/continuous_frechet_test.cpp diff --git a/README.md b/README.md index ea1058f..7ac3be2 100644 --- a/README.md +++ b/README.md @@ -80,6 +80,7 @@ and it then falls back to its default. [engine/configs/geolife.json](engine/conf | `range` | Box size (m), time window (s) and sampling limits for range queries | | `where_at` | Probes per trajectory and max time gap (s) between points to probe | | `where_at_scoring` | Hit radius and tolerance range (m) for scoring where-at queries | +| `frechet` | Fréchet variant for kNN, join and clustering: `Discrete` (kept vertices only, the default) or `Continuous` (polyline distance, unaffected by adding a vertex on an edge) | | `clustering` | Route clustering metric: `enabled`, neighbour threshold `eps_m` (m) and `min_size` per cluster. Reuses the self join when `join.mode` is `Self` with the same eps | | `allocation_alphas` | Allocation exponents to test: budgets follow size^α (`1` = size-proportional, `0` = equal). Exports for α ≠ 1 get an `_a<α>` suffix | diff --git a/engine/CMakeLists.txt b/engine/CMakeLists.txt index 5253712..b0f0f3f 100644 --- a/engine/CMakeLists.txt +++ b/engine/CMakeLists.txt @@ -31,6 +31,7 @@ add_library(engine src/evaluator/knn.cpp src/query/knn.cpp src/query/frechet_helper.cpp + src/query/continuous_frechet.cpp src/evaluator/where_at.cpp src/query/range_query.cpp src/evaluator/range.cpp @@ -72,6 +73,7 @@ add_executable(engine_tests tests/join_test.cpp tests/combined_test.cpp tests/route_clustering_test.cpp + tests/continuous_frechet_test.cpp ) target_link_libraries(engine_tests PRIVATE engine GTest::gtest_main Arrow::arrow_shared Parquet::parquet_shared) gtest_discover_tests(engine_tests) diff --git a/engine/configs/geolife.json b/engine/configs/geolife.json index bf413f9..305a35e 100644 --- a/engine/configs/geolife.json +++ b/engine/configs/geolife.json @@ -33,11 +33,12 @@ }, "join": { "mode": "Self", - "eps_m": 500.0 + "eps_m": 150.0 }, "clustering": { - "enabled": false, + "enabled": true, "eps_m": 150.0, "min_size": 5 - } + }, + "frechet": "Continuous" } diff --git a/engine/include/backend/engine/query/frechet_helper.hpp b/engine/include/backend/engine/query/frechet_helper.hpp index 53e6c03..adf1edb 100644 --- a/engine/include/backend/engine/query/frechet_helper.hpp +++ b/engine/include/backend/engine/query/frechet_helper.hpp @@ -6,6 +6,10 @@ namespace query { +// Discrete compares kept vertices only. Continuous (Alt and Godau) compares the polylines, so +// adding a vertex on an edge never changes it. Both operands of a comparison use the same variant. +enum class FrechetVariant { Discrete, Continuous }; + // Squared instead of std::hypot, which is slow. Fréchet only takes max and min of distances, // Square still has order inline double squared_point_dist(double a_x, double a_y, double b_x, double b_y) { @@ -31,4 +35,15 @@ bool frechet_within(const core::TrajectoryView &a, const core::TrajectoryView &b double endpoint_lower_bound(const core::TrajectoryView &a, const core::TrajectoryView &b); +// Whether the continuous Fréchet distance is at most max_dist, up to a relative 1e-12 +bool continuous_frechet_within(const core::TrajectoryView &a, const core::TrajectoryView &b, + double max_dist); + +// Continuous Fréchet distance within 1 mm, searched in [lower, upper]. Infinity above upper. +double continuous_frechet_distance(const core::TrajectoryView &a, const core::TrajectoryView &b, + double lower, double upper); + +bool frechet_within(const core::TrajectoryView &a, const core::TrajectoryView &b, double max_dist, + FrechetVariant variant); + } // namespace query diff --git a/engine/include/backend/engine/query/join.hpp b/engine/include/backend/engine/query/join.hpp index d5c0921..aec8dd8 100644 --- a/engine/include/backend/engine/query/join.hpp +++ b/engine/include/backend/engine/query/join.hpp @@ -1,6 +1,7 @@ #pragma once #include "backend/engine/core/trajectory_db.hpp" +#include "backend/engine/query/frechet_helper.hpp" #include #include @@ -14,6 +15,7 @@ struct JoinConfig { std::size_t query_count = 10; // 0 = use every trajectory in the query pool double eps_m = 500.0; // join threshold in metres, must be positive JoinMode mode = JoinMode::HeldOut; // Self ignores query_count and the query pool + FrechetVariant frechet = FrechetVariant::Discrete; }; struct JoinQuery { @@ -28,6 +30,7 @@ struct JoinWorkload { std::vector entries; double eps_m = 0.0; JoinMode mode = JoinMode::HeldOut; + FrechetVariant frechet = FrechetVariant::Discrete; // Self entries contain only indices j > i; both operands are evaluated compressed. }; @@ -39,9 +42,11 @@ JoinWorkload generate_join_queries(const core::TrajectoryDb &original, // All unordered database pairs, endpoint-pruned, with reference truth on originals. // No held-out queries or sampling. Nonempty trajectories must retain their endpoints/IDs/order. -JoinWorkload generate_self_join(const core::TrajectoryDb &original, double eps_m = 500.0); +JoinWorkload generate_self_join(const core::TrajectoryDb &original, double eps_m = 500.0, + FrechetVariant variant = FrechetVariant::Discrete); std::vector join(const core::TrajectoryDb &trajectories, - const core::TrajectoryView &query, double eps); + const core::TrajectoryView &query, double eps, + FrechetVariant variant = FrechetVariant::Discrete); } // namespace query diff --git a/engine/include/backend/engine/query/knn.hpp b/engine/include/backend/engine/query/knn.hpp index 87554ad..ba820d0 100644 --- a/engine/include/backend/engine/query/knn.hpp +++ b/engine/include/backend/engine/query/knn.hpp @@ -1,6 +1,7 @@ #pragma once #include "backend/engine/core/trajectory_db.hpp" +#include "backend/engine/query/frechet_helper.hpp" #include #include @@ -17,6 +18,7 @@ struct KnnConfig { std::size_t k = 5; std::size_t query_count = 10; // 0 = use every trajectory in the query pool double p = 0.5; // top-k tau penalty for a pair one list holds neither of, in [0, 1] + FrechetVariant frechet = FrechetVariant::Discrete; }; struct KnnWorkload { @@ -24,6 +26,7 @@ struct KnnWorkload { std::vector entries; std::size_t k; double p = 0.5; + FrechetVariant frechet = FrechetVariant::Discrete; }; // Samples query_count trajectories from query_pool (the held-out split) and stores their k @@ -32,9 +35,11 @@ KnnWorkload generate_knn_queries(const core::TrajectoryDb &original, const core::TrajectoryDb &query_pool, const KnnConfig &config = {}, std::uint64_t seed = 42); -// Indices into trajectories of the k trajectories closest to query by discrete Fréchet -// distance, closest first. Equal distances are ordered by database index. +// Indices into trajectories of the k trajectories closest to query by Fréchet distance, closest +// first. Equal distances are ordered by database index; continuous distances within 1 mm count +// as equal. std::vector knn(const core::TrajectoryDb &trajectories, - const core::TrajectoryView &query, std::size_t k); + const core::TrajectoryView &query, std::size_t k, + FrechetVariant variant = FrechetVariant::Discrete); } // namespace query diff --git a/engine/include/backend/engine/query/route_clustering.hpp b/engine/include/backend/engine/query/route_clustering.hpp index 2161372..f1af806 100644 --- a/engine/include/backend/engine/query/route_clustering.hpp +++ b/engine/include/backend/engine/query/route_clustering.hpp @@ -15,6 +15,7 @@ struct RouteClusteringConfig { bool enabled = false; double eps_m = 500.0; // neighbour threshold, same test as the self join std::size_t min_size = 5; // a cluster needs at least this many trajectories, medoid included + FrechetVariant frechet = FrechetVariant::Discrete; }; using TrajectoryPair = std::pair; @@ -47,10 +48,11 @@ struct RouteClusteringWorkload { std::vector originals; double eps_m = 0.0; std::size_t min_size = 0; + FrechetVariant frechet = FrechetVariant::Discrete; }; -// Clusters the original database. A self join workload over the same trajectories and eps is -// reused instead of recomputing the Fréchet tests on the originals. +// Clusters the original database. A self join workload over the same trajectories, eps and +// Fréchet variant is reused instead of recomputing the Fréchet tests on the originals. RouteClusteringWorkload generate_route_clustering(const core::TrajectoryDb &original, const RouteClusteringConfig &config = {}, const JoinWorkload *self_join = nullptr); diff --git a/engine/src/evaluator/join.cpp b/engine/src/evaluator/join.cpp index 6c1cbd5..b7d7c39 100644 --- a/engine/src/evaluator/join.cpp +++ b/engine/src/evaluator/join.cpp @@ -68,7 +68,8 @@ JoinScores evaluate_join(const core::TrajectoryDb &compressed, // Simplification keeps endpoints, so the original candidates cover every possible match Counts counts; for (const auto j : entry.candidates) { - if (!query::frechet_within(query, compressed.trajectory(j), workload.eps_m)) + if (!query::frechet_within(query, compressed.trajectory(j), workload.eps_m, + workload.frechet)) continue; ++counts.found; if (std::ranges::binary_search(entry.reference, j)) diff --git a/engine/src/evaluator/knn.cpp b/engine/src/evaluator/knn.cpp index 7fbbec6..f4f1a42 100644 --- a/engine/src/evaluator/knn.cpp +++ b/engine/src/evaluator/knn.cpp @@ -102,7 +102,7 @@ KnnEvaluation evaluate_knn(const core::TrajectoryDb &compressed, core::parallel_for(workload.entries.size(), [&](std::size_t i) { const auto &query = workload.entries[i]; const auto trajectory = workload.queries.trajectory(query.query_index); - const auto simplified = query::knn(compressed, trajectory, workload.k); + const auto simplified = query::knn(compressed, trajectory, workload.k, workload.frechet); auto &answer = result.answers[i]; answer.query_trajectory_id = trajectory.trajectory_id; answer.neighbors.reserve(simplified.size()); diff --git a/engine/src/evaluator/route_clustering.cpp b/engine/src/evaluator/route_clustering.cpp index 4a30cd9..3c26bb8 100644 --- a/engine/src/evaluator/route_clustering.cpp +++ b/engine/src/evaluator/route_clustering.cpp @@ -56,7 +56,7 @@ query::RouteClusters compressed_route_clusters(const core::TrajectoryDb &compres core::parallel_for(workload.candidates.size(), [&](std::size_t k) { const auto [i, j] = workload.candidates[k]; matched[k] = query::frechet_within(compressed.trajectory(i), compressed.trajectory(j), - workload.eps_m); + workload.eps_m, workload.frechet); }); std::vector neighbours; for (std::size_t k = 0; k < matched.size(); ++k) { diff --git a/engine/src/pipeline/config.cpp b/engine/src/pipeline/config.cpp index 0785f64..c33ad67 100644 --- a/engine/src/pipeline/config.cpp +++ b/engine/src/pipeline/config.cpp @@ -109,7 +109,7 @@ Config parse_config(const json &data) { validate_keys(data, {"data_path", "export_results", "export_dir", "seed", "budget_shares", "allocation_alphas", "algorithms", "split", "knn", "range", "where_at", - "where_at_scoring", "join", "clustering"}, + "where_at_scoring", "join", "clustering", "frechet"}, "root"); Config config; @@ -144,6 +144,13 @@ Config parse_config(const json &data) { parse_where_at_scoring_config(config, data.value("where_at_scoring", json::object())); parse_join_config(config, data.value("join", json::object())); parse_clustering_config(config, data.value("clustering", json::object())); + // One distance for every Fréchet-based query, so both sides of each comparison agree + if (data.contains("frechet")) { + const auto variant = parse_enum(data["frechet"]); + config.knn.frechet = variant; + config.join.frechet = variant; + config.clustering.frechet = variant; + } return config; } diff --git a/engine/src/query/continuous_frechet.cpp b/engine/src/query/continuous_frechet.cpp new file mode 100644 index 0000000..d73d07a --- /dev/null +++ b/engine/src/query/continuous_frechet.cpp @@ -0,0 +1,200 @@ +#include "backend/engine/query/frechet_helper.hpp" + +#include +#include +#include +#include +#include +#include + +namespace query { + +namespace { + +struct Interval { + double lo = 2.0; + double hi = -1.0; + bool empty() const { return lo > hi; } +}; + +// Parameters t in [0, 1] where |a + t(b - a) - p|^2 <= eps2 +Interval free_interval(double px, double py, double ax, double ay, double bx, double by, + double eps2) { + const double dx = bx - ax, dy = by - ay; + const double wx = ax - px, wy = ay - py; + const double a = dx * dx + dy * dy; + const double start = wx * wx + wy * wy; + const double ex = bx - px, ey = by - py; + const double end = ex * ex + ey * ey; + const bool start_in = start <= eps2; + const bool end_in = end <= eps2; + if (start_in && end_in) + return {0.0, 1.0}; // the disc is convex + if (a == 0.0) + return start_in ? Interval{0.0, 1.0} : Interval{}; + // Closest parameter and squared distance to the line. The cross product avoids the + // cancellation of the quadratic's discriminant on long segments. + const double closest = -(wx * dx + wy * dy) / a; + const double cross = wx * dy - wy * dx; + const double line2 = cross * cross / a; + if (line2 > eps2) + return {}; + const double half = std::sqrt((eps2 - line2) / a); + const double lo = start_in ? 0.0 : std::max(0.0, closest - half); + const double hi = end_in ? 1.0 : std::min(1.0, closest + half); + if (lo > hi) + return {}; + return {lo, hi}; +} + +struct Scratch { + std::vector cur, next; +}; + +// Free-space reachability over the cells of P's segments (columns) against Q's (rows) +bool decide(std::span px, std::span py, std::span qx, + std::span qy, double eps, Scratch &s) { + const std::size_t n = px.size(), m = qx.size(); + if (n == 0 || m == 0) + return false; + const double eps2 = + std::nextafter(eps * eps * (1.0 + 1e-12), std::numeric_limits::infinity()); + const auto d2 = [](double ax, double ay, double bx, double by) { + const double dx = ax - bx, dy = ay - by; + return dx * dx + dy * dy; + }; + if (d2(px[0], py[0], qx[0], qy[0]) > eps2 || + d2(px[n - 1], py[n - 1], qx[m - 1], qy[m - 1]) > eps2) + return false; + // A single point against a polyline: every vertex must be within eps + if (n == 1 || m == 1) { + const auto &ax = n == 1 ? qx : px; + const auto &ay = n == 1 ? qy : py; + const double cx = n == 1 ? px[0] : qx[0]; + const double cy = n == 1 ? py[0] : qy[0]; + for (std::size_t k = 0; k < ax.size(); ++k) + if (d2(ax[k], ay[k], cx, cy) > eps2) + return false; + return true; + } + + const std::size_t segs_q = m - 1; + s.cur.assign(segs_q, Interval{}); + s.next.assign(segs_q, Interval{}); + + // Column 0, left boundary: reachable only as a free prefix from (0, 0) + std::size_t first = segs_q, last = 0; + for (std::size_t j = 0; j < segs_q; ++j) { + const Interval f = free_interval(px[0], py[0], qx[j], qy[j], qx[j + 1], qy[j + 1], eps2); + if (f.empty() || f.lo != 0.0) + break; + s.cur[j] = f; + first = 0; + last = j; + if (f.hi != 1.0) + break; + } + Interval bottom = free_interval(qx[0], qy[0], px[0], py[0], px[1], py[1], eps2); + if (bottom.lo != 0.0) + bottom = {}; + + Interval top_last{}; + for (std::size_t i = 0; i + 1 < n; ++i) { + const double ax = px[i], ay = py[i], bx = px[i + 1], by = py[i + 1]; + std::size_t next_first = segs_q, next_last = 0; + bool any = false; + std::size_t j = bottom.empty() ? first : 0; + if (j >= segs_q) { + if (!bottom.empty()) + j = 0; + else + return false; + } + Interval below = bottom; + for (; j < segs_q; ++j) { + const bool has_left = j >= first && j <= last && !s.cur[j].empty(); + const Interval left = has_left ? s.cur[j] : Interval{}; + if (left.empty() && below.empty()) { + if (j > last) + break; + continue; + } + const double cx = qx[j], cy = qy[j], dx = qx[j + 1], dy = qy[j + 1]; + Interval right = free_interval(bx, by, cx, cy, dx, dy, eps2); + if (below.empty() && !right.empty()) + right.lo = std::max(right.lo, left.lo); + Interval top = free_interval(dx, dy, ax, ay, bx, by, eps2); + if (left.empty() && !top.empty()) + top.lo = std::max(top.lo, below.lo); + if (right.empty()) + right = {}; + if (top.empty()) + top = {}; + s.next[j] = right; + if (!right.empty()) { + if (!any) + next_first = j; + next_last = j; + any = true; + } + below = top; + if (j + 1 == segs_q) + top_last = top; + } + // Bottom boundary for column i + 1 + Interval next_bottom{}; + if (i + 2 < n && !bottom.empty() && bottom.hi == 1.0) { + next_bottom = free_interval(qx[0], qy[0], bx, by, px[i + 2], py[i + 2], eps2); + if (next_bottom.lo != 0.0) + next_bottom = {}; + } + if (i + 2 == n) { + const bool right_corner = + any && next_last == segs_q - 1 && s.next[segs_q - 1].hi == 1.0; + const bool top_corner = j >= segs_q && !top_last.empty() && top_last.hi == 1.0; + return right_corner || top_corner; + } + if (!any && next_bottom.empty()) + return false; + for (std::size_t k = first; k <= last && k < segs_q; ++k) + s.cur[k] = {}; + std::swap(s.cur, s.next); + first = any ? next_first : segs_q; + last = any ? next_last : 0; + bottom = next_bottom; + top_last = {}; + } + return false; +} + +} // namespace + +bool continuous_frechet_within(const core::TrajectoryView &a, const core::TrajectoryView &b, + double max_dist) { + thread_local Scratch scratch; + return decide(a.x, a.y, b.x, b.y, max_dist, scratch); +} + +double continuous_frechet_distance(const core::TrajectoryView &a, const core::TrajectoryView &b, + double lower, double upper) { + constexpr double tolerance = 1e-3; + if (!std::isfinite(upper) || !continuous_frechet_within(a, b, upper)) + return std::numeric_limits::infinity(); + lower = std::max(0.0, std::min(lower, upper)); + while (upper - lower > tolerance) { + const double mid = 0.5 * (lower + upper); + if (continuous_frechet_within(a, b, mid)) + upper = mid; + else + lower = mid; + } + return upper; +} + +bool frechet_within(const core::TrajectoryView &a, const core::TrajectoryView &b, double max_dist, + FrechetVariant variant) { + return variant == FrechetVariant::Continuous ? continuous_frechet_within(a, b, max_dist) + : frechet_within(a, b, max_dist); +} + +} // namespace query diff --git a/engine/src/query/join.cpp b/engine/src/query/join.cpp index 2e378a5..fc648b1 100644 --- a/engine/src/query/join.cpp +++ b/engine/src/query/join.cpp @@ -28,10 +28,11 @@ std::vector endpoint_candidates(const core::TrajectoryDb &trajector } // namespace std::vector join(const core::TrajectoryDb &trajectories, - const core::TrajectoryView &query, double eps) { + const core::TrajectoryView &query, double eps, + FrechetVariant variant) { std::vector result; for (const auto i : endpoint_candidates(trajectories, query, eps)) { - if (frechet_within(query, trajectories.trajectory(i), eps)) + if (frechet_within(query, trajectories.trajectory(i), eps, variant)) result.push_back(i); } return result; @@ -44,7 +45,7 @@ JoinWorkload generate_join_queries(const core::TrajectoryDb &original, throw std::invalid_argument("join: eps_m must be finite and positive"); if (config.mode == JoinMode::Self) - return generate_self_join(original, config.eps_m); + return generate_self_join(original, config.eps_m, config.frechet); std::vector picked(query_pool.size()); std::iota(picked.begin(), picked.end(), 0); @@ -54,7 +55,8 @@ JoinWorkload generate_join_queries(const core::TrajectoryDb &original, picked.resize(config.query_count); } - JoinWorkload results{core::subset(query_pool, picked), {}, config.eps_m}; + JoinWorkload results{ + core::subset(query_pool, picked), {}, config.eps_m, JoinMode::HeldOut, config.frechet}; const std::size_t n = results.queries.size(); if (n == 0 || original.size() == 0) return results; @@ -66,14 +68,15 @@ JoinWorkload generate_join_queries(const core::TrajectoryDb &original, entry.query_index = i; entry.candidates = endpoint_candidates(original, query, results.eps_m); for (const auto j : entry.candidates) { - if (frechet_within(query, original.trajectory(j), results.eps_m)) + if (frechet_within(query, original.trajectory(j), results.eps_m, config.frechet)) entry.reference.push_back(j); } }); return results; } -JoinWorkload generate_self_join(const core::TrajectoryDb &original, double eps_m) { +JoinWorkload generate_self_join(const core::TrajectoryDb &original, double eps_m, + FrechetVariant variant) { if (!std::isfinite(eps_m) || eps_m <= 0.0) throw std::invalid_argument("self join: eps_m must be finite and positive"); core::Simplification endpoints; @@ -85,7 +88,7 @@ JoinWorkload generate_self_join(const core::TrajectoryDb &original, double eps_m n == 1 ? std::vector{0} : std::vector{0, n - 1}; endpoints.add_trajectory(kept); } - JoinWorkload result{core::materialize(original, endpoints), {}, eps_m, JoinMode::Self}; + JoinWorkload result{core::materialize(original, endpoints), {}, eps_m, JoinMode::Self, variant}; result.entries.resize(original.size()); core::parallel_for(original.size(), [&](std::size_t i) { auto &entry = result.entries[i]; @@ -96,7 +99,7 @@ JoinWorkload generate_self_join(const core::TrajectoryDb &original, double eps_m if (endpoint_lower_bound(a, b) > eps_m) continue; entry.candidates.push_back(j); - if (frechet_within(a, b, eps_m)) + if (frechet_within(a, b, eps_m, variant)) entry.reference.push_back(j); } }); diff --git a/engine/src/query/knn.cpp b/engine/src/query/knn.cpp index 8dd4a23..5867486 100644 --- a/engine/src/query/knn.cpp +++ b/engine/src/query/knn.cpp @@ -96,55 +96,66 @@ struct Candidate { std::size_t idx; }; -} // namespace - -struct neighbor { - std::size_t idx; - double squared_frechet_distance; - - bool operator<(const neighbor &other) const { - return std::tie(squared_frechet_distance, idx) < - std::tie(other.squared_frechet_distance, other.idx); - } -}; +double squared_segment_dist(double px, double py, double ax, double ay, double bx, double by) { + // Endpoints are measured directly and interior projections by the cross product, so a point + // on the segment is not pushed off it by reconstructing a + t (b - a) + const double dx = bx - ax, dy = by - ay; + const double length = dx * dx + dy * dy; + const double along = (px - ax) * dx + (py - ay) * dy; + if (length == 0.0 || along <= 0.0) + return squared_point_dist(px, py, ax, ay); + if (along >= length) + return squared_point_dist(px, py, bx, by); + const double cross = (px - ax) * dy - (py - ay) * dx; + return cross * cross / length; +} -std::vector knn(const core::TrajectoryDb &trajectory_db, - const core::TrajectoryView &query_trajectory, std::size_t k) { - if (k == 0) { - throw std::invalid_argument("knn: k = 0, no neighbors"); +// Largest distance from a vertex of a to the polyline b, a lower bound on the continuous +// distance. A vertex stops scanning once it is no farther than the current maximum, and the +// whole scan stops above stop. +double squared_directed_hausdorff(const core::TrajectoryView &a, const core::TrajectoryView &b, + double stop) { + const double inf = std::numeric_limits::infinity(); + if (a.size() == 0 || b.size() == 0) + return inf; + double bound = 0.0; + for (std::size_t i = 0; i < a.size(); ++i) { + double nearest = b.size() == 1 ? squared_point_dist(a.x[i], a.y[i], b.x[0], b.y[0]) : inf; + for (std::size_t j = 0; j + 1 < b.size() && nearest > bound; ++j) + nearest = std::min(nearest, squared_segment_dist(a.x[i], a.y[i], b.x[j], b.y[j], + b.x[j + 1], b.y[j + 1])); + bound = std::max(bound, nearest); + if (bound > stop) + return bound; } + return bound; +} - // 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] = {std::max(squared_endpoint_lower_bound(query_trajectory, reference), - squared_box_extent_lower_bound(query_trajectory, reference)), +// Candidates ordered 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. +class CandidateQueue { + public: + CandidateQueue(const core::TrajectoryDb &trajectory_db, const core::TrajectoryView &query, + std::size_t k) { + m_heap.resize(trajectory_db.size()); + for (std::size_t i = 0; i < trajectory_db.size(); ++i) { + const core::TrajectoryView reference = trajectory_db.trajectory(i); + m_heap[i] = {std::max(squared_endpoint_lower_bound(query, reference), + squared_box_extent_lower_bound(query, 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); - }; - 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()) { + } + std::make_heap(m_heap.begin(), m_heap.end(), farther); + if (k > m_heap.size()) + return; + // The greedy paths bound the k-th distance until exact results tighten it constexpr std::size_t kInitialBoundCandidates = 64; - const auto count = std::min(candidates.size(), std::max(k, kInitialBoundCandidates)); - seed_candidates.reserve(count); + const auto count = std::min(m_heap.size(), std::max(k, kInitialBoundCandidates)); + m_seeds.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(seed_candidates.back().idx)); + m_seeds.push_back(extract()); + const double upper = + squared_greedy_upper_bound(query, trajectory_db.trajectory(m_seeds.back().idx)); if (upper_bounds.size() < k) upper_bounds.push(upper); else if (upper < upper_bounds.top()) { @@ -152,39 +163,68 @@ std::vector knn(const core::TrajectoryDb &trajectory_db, 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()); + m_initial_cutoff = upper_bounds.top(); + } + + bool empty() const { return m_next_seed == m_seeds.size() && m_heap.empty(); } + Candidate next() { return m_next_seed < m_seeds.size() ? m_seeds[m_next_seed++] : extract(); } + // Squared distance the k-th neighbour is known not to exceed, infinity without k candidates + double initial_cutoff() const { return m_initial_cutoff; } + + private: + static bool farther(const Candidate &a, const Candidate &b) { + return std::tie(a.squared_lower_bound, a.idx) > std::tie(b.squared_lower_bound, b.idx); + } + Candidate extract() { + std::pop_heap(m_heap.begin(), m_heap.end(), farther); + const Candidate candidate = m_heap.back(); + m_heap.pop_back(); + return candidate; } + std::vector m_heap; + std::vector m_seeds; + std::size_t m_next_seed = 0; + double m_initial_cutoff = std::numeric_limits::infinity(); +}; + +struct neighbor { + std::size_t idx; + double squared_frechet_distance; + + bool operator<(const neighbor &other) const { + return std::tie(squared_frechet_distance, idx) < + std::tie(other.squared_frechet_distance, other.idx); + } +}; + +std::vector discrete_knn(const core::TrajectoryDb &trajectory_db, + const core::TrajectoryView &query_trajectory, std::size_t k) { + CandidateQueue queue(trajectory_db, query_trajectory, k); + // Include candidates equal to the upper bound, since we still need k results + const double initial_cutoff = + std::nextafter(queue.initial_cutoff(), std::numeric_limits::infinity()); std::priority_queue curr_neighbors; std::vector prev; std::vector curr; - 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(); + while (!queue.empty()) { + const Candidate candidate = queue.next(); const std::size_t i = candidate.idx; - - // 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) { + if (candidate.squared_lower_bound >= curr_max_dist) break; - } core::TrajectoryView reference_trajectory = trajectory_db.trajectory(i); if (std::isfinite(curr_max_dist) && squared_points_to_box_lower_bound(query_trajectory, reference_trajectory, - curr_max_dist) >= curr_max_dist) { + curr_max_dist) >= curr_max_dist) continue; - } - double curr_dist = get_squared_frechet_distance(query_trajectory, reference_trajectory, - curr_max_dist, prev, curr); + const 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) { @@ -204,6 +244,144 @@ std::vector knn(const core::TrajectoryDb &trajectory_db, } return neighbor_ids; } + +// A candidate whose continuous distance is known to lie in [lo, hi] +struct Bracket { + std::size_t idx; + double lo; + double hi; +}; + +// Each distance is narrowed only as far as top-k membership and order need. Distances within +// tolerance of each other count as equal and are ordered by database index, and the same rule +// decides membership at the k-th place. Every prune allows the tolerance as well. +std::vector continuous_knn(const core::TrajectoryDb &trajectory_db, + const core::TrajectoryView &query, std::size_t k) { + constexpr double tolerance = 1e-3; + CandidateQueue queue(trajectory_db, query, k); + std::vector kept; + std::vector scratch; + const auto halve = [&](Bracket &b) { + const double mid = 0.5 * (b.lo + b.hi); + if (continuous_frechet_within(query, trajectory_db.trajectory(b.idx), mid)) + b.hi = mid; + else + b.lo = mid; + }; + const auto kth_upper = [&] { + scratch.clear(); + for (const auto &b : kept) + scratch.push_back(b.hi); + std::nth_element(scratch.begin(), scratch.begin() + static_cast(k - 1), + scratch.end()); + return scratch[k - 1]; + }; + // Narrows overlapping brackets until they are disjoint or within tolerance, then orders runs + // of overlapping brackets, which are ties, by index + const auto rank = [&] { + const auto by_lower = [](const Bracket &b) { return std::pair(b.lo, b.idx); }; + for (bool overlapping = true; overlapping;) { + overlapping = false; + std::ranges::sort(kept, {}, by_lower); + for (std::size_t i = 0; i + 1 < kept.size(); ++i) { + Bracket &a = kept[i], &b = kept[i + 1]; + if (a.hi < b.lo) + continue; + Bracket &wider = a.hi - a.lo >= b.hi - b.lo ? a : b; + if (wider.hi - wider.lo <= tolerance) + continue; + halve(wider); + overlapping = true; + } + } + for (std::size_t start = 0; start < kept.size();) { + std::size_t end = start + 1; + double reach = kept[start].hi; + while (end < kept.size() && kept[end].lo <= reach) + reach = std::max(reach, kept[end++].hi); + std::ranges::sort(kept.begin() + static_cast(start), + kept.begin() + static_cast(end), {}, &Bracket::idx); + start = end; + } + }; + // Drops brackets that are certainly beyond the k-th, refining the widest uncertain one. + // Once only ties remain at the boundary, the ranking decides who stays. + const auto trim = [&] { + while (kept.size() > k) { + const double cutoff = kth_upper(); + std::erase_if(kept, [&](const Bracket &b) { return b.lo > cutoff + tolerance; }); + if (kept.size() <= k) + return; + Bracket *widest = nullptr; + for (auto &b : kept) + if (b.hi >= cutoff && b.hi - b.lo > tolerance && + (!widest || b.hi - b.lo > widest->hi - widest->lo)) + widest = &b; + if (!widest) { + rank(); + kept.resize(k); + return; + } + halve(*widest); + } + }; + + std::vector prev, curr; + while (!queue.empty()) { + const Candidate candidate = queue.next(); + const double cutoff = kept.size() < k ? std::sqrt(queue.initial_cutoff()) : kth_upper(); + const double limit = cutoff + tolerance; + const double squared_limit = limit * limit; + if (candidate.squared_lower_bound > squared_limit) + break; + const auto reference = trajectory_db.trajectory(candidate.idx); + if (squared_points_to_box_lower_bound(query, reference, squared_limit) > squared_limit) + continue; + const double squared_hausdorff = + std::max(squared_directed_hausdorff(query, reference, squared_limit), + squared_directed_hausdorff(reference, query, squared_limit)); + if (squared_hausdorff > squared_limit) + continue; + // The discrete distance bounds the continuous one from above + double hi = limit; + const double squared_discrete = get_squared_frechet_distance( + query, reference, std::nextafter(squared_limit, squared_limit + 1.0), prev, curr); + if (squared_discrete <= squared_limit) + hi = std::min(hi, std::sqrt(squared_discrete)); + else if (!continuous_frechet_within(query, reference, limit)) + continue; + double lo = std::min( + hi, std::max(std::sqrt(candidate.squared_lower_bound), std::sqrt(squared_hausdorff))); + // The distance often equals the Hausdorff bound + if (hi - lo > tolerance) { + if (continuous_frechet_within(query, reference, lo + tolerance)) + hi = lo + tolerance; + else + lo += tolerance; + } + kept.push_back({candidate.idx, lo, hi}); + trim(); + } + + rank(); + std::vector neighbor_ids; + for (const auto &b : kept) + neighbor_ids.push_back(b.idx); + return neighbor_ids; +} + +} // namespace + +std::vector knn(const core::TrajectoryDb &trajectory_db, + const core::TrajectoryView &query_trajectory, std::size_t k, + FrechetVariant variant) { + if (k == 0) { + throw std::invalid_argument("knn: k = 0, no neighbors"); + } + return variant == FrechetVariant::Discrete ? discrete_knn(trajectory_db, query_trajectory, k) + : continuous_knn(trajectory_db, query_trajectory, k); +} + KnnWorkload generate_knn_queries(const core::TrajectoryDb &original, const core::TrajectoryDb &query_pool, const KnnConfig &config, std::uint64_t seed) { @@ -223,11 +401,12 @@ KnnWorkload generate_knn_queries(const core::TrajectoryDb &original, picked.resize(config.query_count); } - KnnWorkload results{core::subset(query_pool, picked), {}, config.k, config.p}; + KnnWorkload results{core::subset(query_pool, picked), {}, config.k, config.p, config.frechet}; results.entries.resize(results.queries.size()); core::parallel_for(results.queries.size(), [&](std::size_t i) { - results.entries[i] = {i, knn(original, results.queries.trajectory(i), config.k)}; + results.entries[i] = { + i, knn(original, results.queries.trajectory(i), config.k, config.frechet)}; }); return results; } diff --git a/engine/src/query/route_clustering.cpp b/engine/src/query/route_clustering.cpp index b173b8f..52d2eeb 100644 --- a/engine/src/query/route_clustering.cpp +++ b/engine/src/query/route_clustering.cpp @@ -87,6 +87,7 @@ RouteClusteringWorkload generate_route_clustering(const core::TrajectoryDb &orig RouteClusteringWorkload workload; workload.eps_m = config.eps_m; workload.min_size = config.min_size; + workload.frechet = config.frechet; workload.originals.reserve(original.size()); for (std::size_t i = 0; i < original.size(); ++i) workload.originals.push_back(endpoints_of(original.trajectory(i))); @@ -105,7 +106,8 @@ RouteClusteringWorkload generate_route_clustering(const core::TrajectoryDb &orig }; std::vector reference; const bool reuse = self_join && self_join->mode == JoinMode::Self && - self_join->eps_m == config.eps_m && same_database(); + self_join->eps_m == config.eps_m && self_join->frechet == config.frechet && + same_database(); if (reuse) { for (const auto &entry : self_join->entries) { for (const auto j : entry.candidates) @@ -114,7 +116,7 @@ RouteClusteringWorkload generate_route_clustering(const core::TrajectoryDb &orig reference.emplace_back(entry.query_index, j); } } else { - const auto join = generate_self_join(original, config.eps_m); + const auto join = generate_self_join(original, config.eps_m, config.frechet); return generate_route_clustering(original, config, &join); } workload.reference = cluster_routes(original, reference, config.min_size); diff --git a/engine/tests/config_test.cpp b/engine/tests/config_test.cpp index c8f31c3..565b6fc 100644 --- a/engine/tests/config_test.cpp +++ b/engine/tests/config_test.cpp @@ -207,3 +207,17 @@ TEST(Config, RejectsInvalidClustering) { std::string(R"({"data_path": "d", "clustering": {"min_size": )") + value + "}}")), std::invalid_argument); } + +TEST(Config, FrechetVariantAppliesToEveryFrechetQuery) { + const auto config = pipeline::load_config( + write_config("frechet.json", R"({"data_path": "d", "frechet": "Continuous"})")); + EXPECT_EQ(config.knn.frechet, query::FrechetVariant::Continuous); + EXPECT_EQ(config.join.frechet, query::FrechetVariant::Continuous); + EXPECT_EQ(config.clustering.frechet, query::FrechetVariant::Continuous); + const auto defaults = + pipeline::load_config(write_config("frechet_default.json", R"({"data_path": "d"})")); + EXPECT_EQ(defaults.join.frechet, query::FrechetVariant::Discrete); + EXPECT_THROW(pipeline::load_config( + write_config("frechet_bad.json", R"({"data_path": "d", "frechet": "Weak"})")), + std::invalid_argument); +} diff --git a/engine/tests/continuous_frechet_test.cpp b/engine/tests/continuous_frechet_test.cpp new file mode 100644 index 0000000..3aee28c --- /dev/null +++ b/engine/tests/continuous_frechet_test.cpp @@ -0,0 +1,307 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "backend/engine/core/projection.hpp" +#include "backend/engine/core/trajectory_db.hpp" +#include "backend/engine/evaluator/join.hpp" +#include "backend/engine/evaluator/route_clustering.hpp" +#include "backend/engine/query/frechet_helper.hpp" +#include "backend/engine/query/join.hpp" +#include "backend/engine/query/knn.hpp" +#include "backend/engine/query/route_clustering.hpp" + +namespace { + +using Line = std::vector>; // (x, y) in local metres + +core::TrajectoryDb make_db(const std::vector &lines) { + std::vector offsets{0}; + std::vector ids; + std::vector t; + std::vector lat; + std::vector lon; + for (std::size_t i = 0; i < lines.size(); ++i) { + ids.push_back(static_cast(i)); + for (std::size_t k = 0; k < lines[i].size(); ++k) { + const auto ll = core::to_lat_lon(lines[i][k].first, lines[i][k].second); + t.push_back(static_cast(k)); + lat.push_back(ll.lat); + lon.push_back(ll.lon); + } + offsets.push_back(t.size()); + } + return {std::move(offsets), std::move(ids), std::move(t), std::move(lat), std::move(lon)}; +} + +double continuous(const core::TrajectoryDb &db, std::size_t a, std::size_t b) { + return query::continuous_frechet_distance(db.trajectory(a), db.trajectory(b), 0.0, 1e7); +} + +std::vector random_lines(std::mt19937_64 &rng, std::size_t count, std::size_t max_len) { + std::uniform_real_distribution start(-300.0, 300.0); + std::normal_distribution step(40.0, 60.0); + std::uniform_int_distribution length(2, max_len); + std::vector lines(count); + for (auto &line : lines) { + double x = start(rng), y = start(rng); + for (std::size_t k = 0, n = length(rng); k < n; ++k) { + line.emplace_back(x, y); + x += step(rng); + y += step(rng) * 0.5; + } + } + return lines; +} + +} // namespace + +TEST(ContinuousFrechet, AddingAVertexOnAnEdgeChangesNothing) { + // Discrete Fréchet jumps to half the segment once one side gains its true midpoint + const auto db = make_db({{{0, 0}, {1000, 0}}, {{0, 0}, {500, 0}, {1000, 0}}}); + EXPECT_NEAR(continuous(db, 0, 1), 0.0, 1e-2); + EXPECT_NEAR(query::frechet_distance(db.trajectory(0), db.trajectory(1)), 500.0, 1e-2); + EXPECT_TRUE(query::continuous_frechet_within(db.trajectory(0), db.trajectory(1), 1.0)); + EXPECT_FALSE(query::frechet_within(db.trajectory(0), db.trajectory(1), 1.0)); +} + +TEST(ContinuousFrechet, KnownDistances) { + // A tent of height 1 against its base, parallel segments 3 apart, and a backtrack of 0.5 + const auto db = make_db({{{0, 0}, {1, 1}, {2, 0}}, + {{0, 0}, {2, 0}}, + {{0, 3}, {2, 3}}, + {{0, 0}, {1.5, 0}, {1, 0}, {2, 0}}}); + EXPECT_NEAR(continuous(db, 0, 1), 1.0, 1e-2); + EXPECT_NEAR(continuous(db, 1, 2), 3.0, 1e-2); + EXPECT_NEAR(continuous(db, 3, 1), 0.25, 1e-2); +} + +TEST(ContinuousFrechet, NeverAboveDiscreteAndInvariantUnderSubdivision) { + std::mt19937_64 rng(5); + const auto lines = random_lines(rng, 120, 12); + std::vector subdivided; + for (const auto &line : lines) { + Line s{line.front()}; + for (std::size_t k = 1; k < line.size(); ++k) { + const auto [ax, ay] = line[k - 1]; + const auto [bx, by] = line[k]; + s.emplace_back(0.5 * (ax + bx), 0.5 * (ay + by)); + s.push_back(line[k]); + } + subdivided.push_back(s); + } + std::vector both = lines; + both.insert(both.end(), subdivided.begin(), subdivided.end()); + const auto db = make_db(both); + const std::size_t n = lines.size(); + for (std::size_t a = 0; a + 1 < n; a += 2) { + const std::size_t b = a + 1; + const double c = continuous(db, a, b); + EXPECT_LE(c, query::frechet_distance(db.trajectory(a), db.trajectory(b)) + 1e-2); + EXPECT_NEAR(c, continuous(db, a + n, b + n), 2e-2); + EXPECT_NEAR(c, continuous(db, b, a), 2e-2); + } +} + +TEST(ContinuousFrechet, DecisionAgreesWithDistance) { + std::mt19937_64 rng(9); + const auto db = make_db(random_lines(rng, 60, 10)); + for (std::size_t a = 0; a + 1 < db.size(); ++a) { + const double d = continuous(db, a, a + 1); + EXPECT_TRUE( + query::continuous_frechet_within(db.trajectory(a), db.trajectory(a + 1), d + 0.01)); + if (d > 0.02) + EXPECT_FALSE( + query::continuous_frechet_within(db.trajectory(a), db.trajectory(a + 1), d - 0.01)); + } +} + +TEST(ContinuousFrechet, KnnMatchesBruteForce) { + std::mt19937_64 rng(13); + const auto db = make_db(random_lines(rng, 80, 10)); + const auto queries = make_db(random_lines(rng, 12, 10)); + for (std::size_t q = 0; q < queries.size(); ++q) { + std::vector> all; + for (std::size_t i = 0; i < db.size(); ++i) + all.emplace_back(query::continuous_frechet_distance(queries.trajectory(q), + db.trajectory(i), 0.0, 1e7), + i); + std::sort(all.begin(), all.end()); + const auto got = + query::knn(db, queries.trajectory(q), 5, query::FrechetVariant::Continuous); + ASSERT_EQ(got.size(), 5u); + // Distances are exact to 1 mm, so compare the distance of each rank + for (std::size_t r = 0; r < 5; ++r) { + const double d = query::continuous_frechet_distance(queries.trajectory(q), + db.trajectory(got[r]), 0.0, 1e7); + EXPECT_NEAR(d, all[r].first, 5e-3); + } + } +} + +TEST(ContinuousFrechet, KnnOrdersEqualDistancesByIndex) { + // Copies of one line at the same distance, plus a nearer and a farther line + const Line line{{0, 10}, {500, 10}, {1000, 10}}; + const auto db = make_db({line, {{0, 50}, {1000, 50}}, line, {{0, 2}, {1000, 2}}, line}); + const auto query = make_db({{{0, 0}, {1000, 0}}}); + const auto got = query::knn(db, query.trajectory(0), 4, query::FrechetVariant::Continuous); + EXPECT_EQ(got, (std::vector{3, 0, 2, 4})); +} + +TEST(ContinuousFrechet, KnnMatchesBruteForceOnDenseTrajectories) { + // Long, densely sampled lines take the discrete upper-bound path + std::mt19937_64 rng(17); + const auto db = make_db(random_lines(rng, 60, 80)); + const auto queries = make_db(random_lines(rng, 8, 80)); + for (std::size_t q = 0; q < queries.size(); ++q) { + std::vector all; + for (std::size_t i = 0; i < db.size(); ++i) + all.push_back(query::continuous_frechet_distance(queries.trajectory(q), + db.trajectory(i), 0.0, 1e7)); + std::ranges::sort(all); + const auto got = + query::knn(db, queries.trajectory(q), 7, query::FrechetVariant::Continuous); + ASSERT_EQ(got.size(), 7u); + for (std::size_t r = 0; r < got.size(); ++r) + EXPECT_NEAR(query::continuous_frechet_distance(queries.trajectory(q), + db.trajectory(got[r]), 0.0, 1e7), + all[r], 5e-3); + } +} + +TEST(ContinuousFrechet, SelfJoinMatchesBruteForceAndScoresOriginalPerfectly) { + std::mt19937_64 rng(21); + const auto db = make_db(random_lines(rng, 60, 8)); + const double eps = 150.0; + const auto workload = query::generate_self_join(db, eps, query::FrechetVariant::Continuous); + std::size_t expected = 0, got = 0; + for (std::size_t i = 0; i < db.size(); ++i) { + got += workload.entries[i].reference.size(); + for (std::size_t j = i + 1; j < db.size(); ++j) + expected += query::continuous_frechet_within(db.trajectory(i), db.trajectory(j), eps); + } + EXPECT_EQ(got, expected); + const auto scores = evaluator::evaluate_join(db, workload); + if (expected > 0) + EXPECT_DOUBLE_EQ(scores.f1, 1.0); +} + +TEST(ContinuousFrechet, ClusteringReusesOnlyASelfJoinOfTheSameVariant) { + // Three copies of a line and a fourth line that matches them only under continuous distance + const Line line{{0, 0}, {1000, 0}}; + const Line subdivided{{0, 0}, {500, 0}, {1000, 0}}; + const auto db = make_db({line, line, line, subdivided}); + const auto discrete_join = query::generate_self_join(db, 100.0); + const query::RouteClusteringConfig config{.enabled = true, + .eps_m = 100.0, + .min_size = 4, + .frechet = query::FrechetVariant::Continuous}; + const auto workload = query::generate_route_clustering(db, config, &discrete_join); + ASSERT_EQ(workload.reference.members.size(), 1u); + EXPECT_EQ(workload.reference.members[0].size(), 4u); + EXPECT_DOUBLE_EQ(evaluator::evaluate_route_clustering(db, workload).score, 1.0); +} + +TEST(ContinuousFrechet, KnnFindsAnExactMatchFarFromTheOrigin) { + const auto db = make_db({{{-10000.123, 12345.678}, {0.01, 0.02}}}); + EXPECT_EQ(query::knn(db, db.trajectory(0), 1, query::FrechetVariant::Continuous), + (std::vector{0})); +} + +TEST(ContinuousFrechet, KnnTiesFromDifferentGeometriesAreConsistentAcrossK) { + // A tent of height 10 and a parallel line 10 m away are equally far from the query + const Line tent{{0, 0}, {500, 10}, {1000, 0}}; + const Line parallel{{0, 10}, {1000, 10}}; + const auto query = make_db({{{0, 0}, {1000, 0}}}); + for (const auto &[lines, first] : {std::pair{std::vector{tent, parallel}, 0u}, + std::pair{std::vector{parallel, tent}, 0u}}) { + const auto db = make_db(lines); + const auto k = query::FrechetVariant::Continuous; + EXPECT_EQ(query::knn(db, query.trajectory(0), 1, k), (std::vector{first})); + EXPECT_EQ(query::knn(db, query.trajectory(0), 2, k), (std::vector{0, 1})); + } +} + +// A view straight from local coordinates, without the lat/lon round trip +struct RawLine { + std::vector t; + std::vector x, y; + explicit RawLine(const Line &line) { + for (std::size_t k = 0; k < line.size(); ++k) { + t.push_back(static_cast(k)); + x.push_back(line[k].first); + y.push_back(line[k].second); + } + } + core::TrajectoryView view() const { return {0, t, {}, {}, x, y}; } +}; + +TEST(ContinuousFrechet, LongSegmentWithAVertexOnItIsAtDistanceZeroExactCoordinates) { + const double bx = -81677.87045, by = -85151.97823; + for (const double f : {0.75641311848, 0.5, 0.1, 0.999, 0.333333333}) { + const RawLine a({{0, 0}, {bx, by}}); + const RawLine b({{0, 0}, {f * bx, f * by}, {bx, by}}); + EXPECT_LT(query::continuous_frechet_distance(a.view(), b.view(), 0.0, 1e6), 1e-3) << f; + EXPECT_TRUE(query::continuous_frechet_within(a.view(), b.view(), 1e-3)) << f; + EXPECT_TRUE(query::continuous_frechet_within(b.view(), a.view(), 1e-3)) << f; + } +} + +TEST(ContinuousFrechet, LongSegmentWithAVertexOnItIsAtDistanceZero) { + // Built in projected coordinates, so the inserted vertex lies on the segment + const auto ends = make_db({{{0, 0}, {-81677.87045, -85151.97823}}}).trajectory(0); + const double ax = ends.x[0], ay = ends.y[0], bx = ends.x[1], by = ends.y[1]; + for (const double f : {0.75641311848, 0.5, 0.1, 0.999}) { + const auto db = make_db( + {{{ax, ay}, {bx, by}}, {{ax, ay}, {ax + f * (bx - ax), ay + f * (by - ay)}, {bx, by}}}); + EXPECT_LT(continuous(db, 0, 1), 1e-3) << f; + EXPECT_TRUE(query::continuous_frechet_within(db.trajectory(0), db.trajectory(1), 1e-3)) + << f; + } +} + +TEST(ContinuousFrechet, NearTangentSegmentIsJustWithinItsDistance) { + // The vertex (500, d) touches the disc of radius d around the segment's closest point + const auto db = make_db({{{0, 0}, {1000, 0}}, {{0, 0}, {500, 37.5}, {1000, 0}}}); + EXPECT_NEAR(continuous(db, 0, 1), 37.5, 2e-3); + EXPECT_TRUE(query::continuous_frechet_within(db.trajectory(0), db.trajectory(1), 37.5 + 1e-3)); + EXPECT_FALSE(query::continuous_frechet_within(db.trajectory(0), db.trajectory(1), 37.5 - 1e-2)); +} + +TEST(ContinuousFrechet, AgreesWithDiscreteDistanceOnDenselySubdividedLines) { + // An independent check: on lines subdivided into steps of at most s, the discrete distance + // is at least the continuous one and at most s above it + std::mt19937_64 rng(31); + const auto lines = random_lines(rng, 40, 6); + const double step = 2.0; + std::vector dense; + for (const auto &line : lines) { + Line d{line.front()}; + for (std::size_t k = 1; k < line.size(); ++k) { + const auto [ax, ay] = line[k - 1]; + const auto [bx, by] = line[k]; + const auto parts = + static_cast(std::ceil(std::hypot(bx - ax, by - ay) / step)); + for (std::size_t s = 1; s <= parts; ++s) { + const double f = static_cast(s) / static_cast(parts); + d.emplace_back(ax + f * (bx - ax), ay + f * (by - ay)); + } + } + dense.push_back(d); + } + const auto sparse_db = make_db(lines); + const auto dense_db = make_db(dense); + for (std::size_t a = 0; a + 1 < lines.size(); a += 2) { + const double c = continuous(sparse_db, a, a + 1); + const double d = + query::frechet_distance(dense_db.trajectory(a), dense_db.trajectory(a + 1)); + EXPECT_LE(c, d + 2e-3); + EXPECT_GE(c, d - step - 1e-2); + } +} From 74edf3eb5025a7f5e7e5107c02ee659c82c2e0a6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristoffer=20Vesterg=C3=A5rd=20Johansson?= <108988684+Kristofferjoha@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:11:06 +0100 Subject: [PATCH 3/3] Update geolife.json --- engine/configs/geolife.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engine/configs/geolife.json b/engine/configs/geolife.json index 305a35e..55d7423 100644 --- a/engine/configs/geolife.json +++ b/engine/configs/geolife.json @@ -1,5 +1,5 @@ { - "data_path": "data/geolife_gap60.parquet", + "data_path": "data/geolife.parquet", "export_results": true, "export_dir": "data/output/geolife", "seed": 42,