Cuda greedy resolver

README.md

@@ -13,7 +13,7 @@ Demonstrator tracking chain for accelerators.
 |                    | Seed finding           | ✅  | ✅   | ✅   | ✅     | ⚪     | ⚪      |
 |                    | Track param estimation | ✅  | ✅   | ✅   | ✅     | ⚪     | ⚪      |
 | **Track finding**  | Combinatorial KF       | ✅  | ✅   | ✅   | 🟡     | ⚪     | ⚪      |
-| **Ambiguity resolution**  | Greedy resolver   | ✅  | ⚪   |  ⚪  | ⚪     | ⚪     | ⚪      |
+| **Ambiguity resolution**  | Greedy resolver   | ✅  | ✅   |  ⚪  | ⚪     | ⚪     | ⚪      |
 | **Track fitting**  | KF                     | ✅  | ✅   | 🟡   | ⚪     | ⚪     | ⚪      |
 
 ✅: exists, 🟡: work started, ⚪: work not started yet

core/include/traccc/ambiguity_resolution/legacy/greedy_ambiguity_resolution_algorithm.hpp

@@ -73,8 +73,8 @@ class greedy_ambiguity_resolution_algorithm
         ///
         /// There is no (track_id) in this algorithm, only (track_index).
 
-        /// Associates each track_index with the track's chi2 value
-        std::vector<traccc::scalar> track_chi2;
+        /// Associates each track_index with the track's p-value
+        std::vector<traccc::scalar> track_pval;
 
         /// Associates each track_index to the track's (measurement_id)s list
         std::vector<std::vector<std::size_t>> measurements_per_track;

core/src/ambiguity_resolution/greedy_ambiguity_resolution_algorithm.cpp

@@ -29,6 +29,8 @@ greedy_ambiguity_resolution_algorithm::operator()(
 
     const std::size_t n_tracks = track_candidates.size();
 
+    // Make sure there is no duplicate measurement
+
     // Make the output container
     track_candidate_container_types::host output{&m_mr.main};
 
@@ -43,14 +45,14 @@ greedy_ambiguity_resolution_algorithm::operator()(
     std::vector<unsigned int> accepted_ids(n_tracks);
     std::iota(accepted_ids.begin(), accepted_ids.end(), 0);
 
-    // Make measurement ID, chi2 and n_measurement vector
+    // Make measurement ID, pval and n_measurement vector
     std::vector<std::vector<std::size_t>> meas_ids(n_tracks);
-    std::vector<traccc::scalar> chi_squares(n_tracks);
-    std::vector<std::size_t> n_meas(n_tracks);
+    std::vector<traccc::scalar> pvals(n_tracks);
+    std::vector<std::size_t> n_meas(n_tracks, 0);
 
     for (unsigned int i = 0; i < n_tracks; i++) {
-        // Fill the chi-squares vectors
-        chi_squares[i] = track_candidates.at(i).header.trk_quality.chi2;
+        // Fill the pval vectors
+        pvals[i] = track_candidates.at(i).header.trk_quality.pval;
 
         const auto& candidates = track_candidates.at(i).items;
         const unsigned int n_cands = candidates.size();
@@ -90,7 +92,19 @@ greedy_ambiguity_resolution_algorithm::operator()(
     tracks_per_measurement.resize(unique_meas.size());
 
     for (const auto& i : accepted_ids) {
-        for (const auto& meas_id : meas_ids[i]) {
+
+        // Make sure there is no duplicate
+        /*
+        assert(meas_ids[i].size() != std::unordered_set<std::size_t>(
+                                         meas_ids[i].begin(), meas_ids[i].end())
+                                         .size());
+        */
+        std::unordered_set<std::size_t> deduplicated_ids(meas_ids[i].begin(),
+                                                         meas_ids[i].end());
+
+        for (const auto& meas_id : deduplicated_ids) {
+
+            // for (const auto& meas_id : meas_ids[i]) {
             const auto it = std::lower_bound(unique_meas.begin(),
                                              unique_meas.end(), meas_id);
             assert(it != unique_meas.end());
@@ -101,14 +115,15 @@ greedy_ambiguity_resolution_algorithm::operator()(
     }
 
     // Count the number of shared measurements
-    std::vector<unsigned int> n_shared(n_tracks);
+    std::vector<unsigned int> n_shared(n_tracks, 0);
     for (const auto& i : accepted_ids) {
         for (const auto& meas_id : meas_ids[i]) {
             const auto it = std::lower_bound(unique_meas.begin(),
                                              unique_meas.end(), meas_id);
             assert(it != unique_meas.end());
             const std::size_t unique_meas_idx = static_cast<std::size_t>(
                 std::distance(unique_meas.begin(), it));
+
             if (tracks_per_measurement[unique_meas_idx].size() > 1) {
                 n_shared[i]++;
             }
@@ -122,32 +137,36 @@ greedy_ambiguity_resolution_algorithm::operator()(
                         static_cast<traccc::scalar>(n_meas[i]);
     }
 
-    // Sort the track id with rel_shared and chi2 to find the worst track fast
+    // Sort the track id with rel_shared and pval to find the worst track
+    // fast
     std::vector<unsigned int> sorted_ids = accepted_ids;
 
-    auto track_comparator = [&rel_shared, &chi_squares](unsigned int a,
-                                                        unsigned int b) {
+    auto track_comparator = [&rel_shared, &pvals](unsigned int a,
+                                                  unsigned int b) {
         if (rel_shared[a] != rel_shared[b]) {
             return rel_shared[a] < rel_shared[b];
         }
-        return chi_squares[a] < chi_squares[b];
+        return pvals[a] > pvals[b];
     };
     std::sort(sorted_ids.begin(), sorted_ids.end(), track_comparator);
 
     // Iterate over tracks
     for (unsigned int iter = 0; iter < m_config.max_iterations; iter++) {
+
         // Terminate if there are no tracks to iterate
         if (accepted_ids.empty()) {
             break;
         }
 
-        unsigned int max_shared{0u};
+        unsigned int max_shared(0);
         for (const auto& i : accepted_ids) {
-            if (n_shared[i] > max_shared)
+            if (n_shared[i] > max_shared) {
                 max_shared = n_shared[i];
+            }
         }
 
-        // Terminate if the max shared measurements is less than the cut value
+        // Terminate if the max shared measurements is less than the cut
+        // value
         if (max_shared < m_config.max_shared_meas) {
             break;
         }
@@ -164,7 +183,14 @@ greedy_ambiguity_resolution_algorithm::operator()(
         // Pop the worst (rejected) id from the sorted ids
         sorted_ids.pop_back();
 
-        const auto& meas_ids_to_remove = meas_ids[worst_track];
+        std::unordered_set<std::size_t> seen;
+        std::vector<std::size_t> meas_ids_to_remove;
+        for (const auto& id : meas_ids[worst_track]) {
+            if (seen.insert(id).second) {
+                meas_ids_to_remove.push_back(id);
+            }
+        }
+
         for (const auto& id : meas_ids_to_remove) {
             const auto it =
                 std::lower_bound(unique_meas.begin(), unique_meas.end(), id);
@@ -174,29 +200,39 @@ greedy_ambiguity_resolution_algorithm::operator()(
 
             auto& tracks = tracks_per_measurement[unique_meas_idx];
 
-            // Remove the worst (rejected) id from the tracks associated with
-            // measurement
+            if (tracks.empty()) {
+                continue;
+            }
+
+            // Remove the worst (rejected) id from the tracks associated
+            // with measurement
             const auto it2 =
                 std::lower_bound(tracks.begin(), tracks.end(), worst_track);
             assert(it2 != tracks.end() && *it2 == worst_track);
+
             tracks.erase(it2);
 
             // If there is only one track associated with measurement, the
             // number of shared measurement can be reduced by one
             if (tracks.size() == 1) {
                 const auto tid = tracks[0];
-                n_shared[tid]--;
+
+                n_shared[tid] -= static_cast<unsigned int>(
+                    std::count(meas_ids[tid].begin(), meas_ids[tid].end(), id));
                 rel_shared[tid] = static_cast<traccc::scalar>(n_shared[tid]) /
                                   static_cast<traccc::scalar>(n_meas[tid]);
-                // Reposition the track to the next of the worse track which is
-                // firstly found during the reverse iteration
+                // Reposition the track to the next of the worse track which
+                // is firstly found during the reverse iteration
                 const auto it3 =
                     std::find(sorted_ids.begin(), sorted_ids.end(), tid);
                 assert(it3 != sorted_ids.end());
                 const auto it4 = std::lower_bound(sorted_ids.begin(), it3, tid,
                                                   track_comparator);
-                sorted_ids.erase(it3);
-                sorted_ids.insert(it4, tid);
+
+                if (it3 != it4) {
+                    sorted_ids.erase(it3);
+                    sorted_ids.insert(it4, tid);
+                }
             }
         }
 

core/src/ambiguity_resolution/legacy/greedy_ambiguity_resolution_algorithm.cpp

@@ -115,8 +115,8 @@ void greedy_ambiguity_resolution_algorithm::compute_initial_state(
             }
         }
 
-        // Add this track chi2 value
-        state.track_chi2.push_back(fit_res.trk_quality.chi2);
+        // Add this track pval value
+        state.track_pval.push_back(fit_res.trk_quality.pval);
         // Add all the (measurement_id)s of this track
         state.measurements_per_track.push_back(std::move(measurements));
         // Initially, every track is in the selected_track list. They will later
@@ -177,7 +177,7 @@ void greedy_ambiguity_resolution_algorithm::resolve(state_t& state) const {
 
     /// Compares two tracks in order to find the one which should be evicted.
     /// First we compare the relative amount of shared measurements. If that is
-    /// indecisive we use the chi2.
+    /// indecisive we use the pval.
     auto track_comperator = [&state](std::pair<std::size_t, std::size_t> a,
                                      std::pair<std::size_t, std::size_t> b) {
         /// Helper to calculate the relative amount of shared measurements.
@@ -191,7 +191,7 @@ void greedy_ambiguity_resolution_algorithm::resolve(state_t& state) const {
             return relative_shared_measurements(a.second) <
                    relative_shared_measurements(b.second);
         }
-        return state.track_chi2[a.second] < state.track_chi2[b.second];
+        return state.track_pval[a.second] > state.track_pval[b.second];
     };
 
     std::size_t iteration_count = 0;

device/common/include/traccc/ambiguity_resolution/device/count_shared_measurements.hpp

@@ -0,0 +1,67 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2025 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Local include(s).
+#include "traccc/device/global_index.hpp"
+
+// Project include(s)
+#include "traccc/definitions/qualifiers.hpp"
+
+// VecMem include(s).
+#include <vecmem/containers/data/jagged_vector_view.hpp>
+#include <vecmem/containers/data/vector_view.hpp>
+#include <vecmem/containers/device_vector.hpp>
+#include <vecmem/containers/jagged_device_vector.hpp>
+
+namespace traccc::device {
+
+/// (Event Data) Payload for the @c traccc::device::count_shared_measurements
+/// function
+struct count_shared_measurements_payload {
+
+    /**
+     * @brief View object to the accepted ids
+     */
+    vecmem::data::vector_view<const unsigned int> accepted_ids_view;
+
+    /**
+     * @brief View object to the vector of measured ids per track
+     */
+    vecmem::data::jagged_vector_view<const std::size_t> meas_ids_view;
+
+    /**
+     * @brief View object to the unique measurement ids
+     */
+    vecmem::data::vector_view<const std::size_t> unique_meas_view;
+
+    /**
+     * @brief View object to the tracks per measurement
+     */
+    vecmem::data::vector_view<const unsigned int>
+        n_accepted_tracks_per_measurement_view;
+
+    /**
+     * @brief View object to the number of shared measurements
+     */
+    vecmem::data::vector_view<unsigned int> n_shared_view;
+};
+
+/// Function used for fill_vectors
+///
+/// @param[in] globalIndex   The index of the current thread
+/// @param[inout] payload      The function call payload
+///
+TRACCC_HOST_DEVICE inline void count_shared_measurements(
+    global_index_t globalIndex,
+    const count_shared_measurements_payload& payload);
+
+}  // namespace traccc::device
+
+// Include the implementation.
+#include "./impl/count_shared_measurements.ipp"

device/common/include/traccc/ambiguity_resolution/device/fill_track_candidates.hpp

@@ -0,0 +1,61 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2025 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Local include(s).
+#include "traccc/device/global_index.hpp"
+
+// Project include(s)
+#include "traccc/definitions/qualifiers.hpp"
+#include "traccc/edm/track_candidate.hpp"
+
+// VecMem include(s).
+#include <vecmem/containers/data/jagged_vector_view.hpp>
+#include <vecmem/containers/data/vector_view.hpp>
+#include <vecmem/containers/device_vector.hpp>
+#include <vecmem/containers/jagged_device_vector.hpp>
+
+namespace traccc::device {
+
+/// (Event Data) Payload for the @c traccc::device::fill_track_candidates
+/// function
+struct fill_track_candidates_payload {
+
+    /**
+     * @brief View object to the input track candidates
+     */
+    track_candidate_container_types::const_view track_candidates_view;
+
+    /**
+     * @brief The number of accepted tracks
+     */
+    unsigned int n_accepted;
+
+    /**
+     * @brief View object to the sorted ids
+     */
+    vecmem::data::vector_view<const unsigned int> sorted_ids_view;
+
+    /**
+     * @brief View object to the output track candidates
+     */
+    track_candidate_container_types::view res_track_candidates_view;
+};
+
+/// Function used for fill_vectors
+///
+/// @param[in] globalIndex   The index of the current thread
+/// @param[inout] payload      The function call payload
+///
+TRACCC_HOST_DEVICE inline void fill_track_candidates(
+    global_index_t globalIndex, const fill_track_candidates_payload& payload);
+
+}  // namespace traccc::device
+
+// Include the implementation.
+#include "./impl/fill_track_candidates.ipp"