clustering

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

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)

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
     }
 }

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;
 };

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 <limits>
 #include <span>
@@ -19,6 +21,7 @@ struct Workloads {
     std::span<const query::WhereAtQuery> where_at = {};
     std::span<const query::RangeQuery> 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<double>::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;

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 <cstddef>
+#include <limits>
+
+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<double>::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

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.

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 <cstddef>
+#include <cstdint>
+#include <span>
+#include <utility>
+#include <vector>
+
+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<std::size_t, std::size_t>;
+
+// 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<std::ptrdiff_t> label;             // cluster per trajectory, -1 = unclustered
+    std::vector<std::vector<std::size_t>> members; // members[c][0] is the medoid
+};
+
+RouteClusters cluster_routes(const core::TrajectoryDb &trajectories,
+                             std::span<const TrajectoryPair> 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<TrajectoryPair> candidates; // endpoint-pruned pairs, i < j
+    RouteClusters reference;
+    std::vector<TrajectoryEndpoints> 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

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;
 }
 

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

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 <algorithm>
+#include <cstddef>
+#include <stdexcept>
+#include <unordered_map>
+#include <vector>
+
+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<std::ptrdiff_t, std::size_t> 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<double>(shared) / static_cast<double>(size));
+        }
+        sum += best;
+    }
+    return sum / static_cast<double>(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<char> 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<query::TrajectoryPair> 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

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<std::size_t>() < 2)
+            throw std::invalid_argument("config: clustering.min_size must be an integer >= 2");
+        config.clustering.min_size = min_size.get<std::size_t>();
+    }
+}
+
 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;
 }