Add sanity checks for contiguity and orderedness

core/include/traccc/finding/finding_algorithm.ipp

@@ -7,6 +7,8 @@
 
 // Project include(s).
 #include "traccc/finding/candidate_link.hpp"
+#include "traccc/sanity/contiguous_on.hpp"
+#include "traccc/utils/projections.hpp"
 
 // detray include(s).
 #include "detray/geometry/barcode.hpp"
@@ -29,6 +31,10 @@ finding_algorithm<stepper_t, navigator_t>::operator()(
      * Measurement Operations
      *****************************************************************/
 
+    // Check contiguity of the measurements
+    assert(host::is_contiguous_on(measurement_module_projection(),
+                                  vecmem::get_data(measurements)));
+
     // Get copy of barcode uniques
     std::vector<measurement> uniques;
     uniques.resize(measurements.size());

core/include/traccc/sanity/contiguous_on.hpp

@@ -0,0 +1,94 @@
+/**
+ * traccc library, part of the ACTS project (R&D line)
+ *
+ * (c) 2024 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s).
+#include "traccc/definitions/concepts.hpp"
+
+// VecMem include(s).
+#include <memory>
+#include <vecmem/containers/data/vector_view.hpp>
+#include <vecmem/containers/device_vector.hpp>
+#include <vecmem/memory/memory_resource.hpp>
+#include <vecmem/utils/copy.hpp>
+
+// System include
+#if __cpp_concepts >= 201907L
+#include <concepts>
+#endif
+#include <memory>
+#include <unordered_set>
+
+namespace traccc::host {
+/**
+ * @brief Sanity check that a given vector is contiguous on a given projection.
+ *
+ * For a vector $v$ to be contiguous on a projection $\pi$, it must be the case
+ * that for all indices $i$ and $j$, if $v_i = v_j$, then all indices $k$
+ * between $i$ and $j$, $v_i = v_j = v_k$.
+ *
+ * @note This function runs in O(n^2) time.
+ *
+ * @tparam P The type of projection $\pi$, a callable which returns some
+ * comparable type.
+ * @tparam T The type of the vector.
+ * @param projection A projection object of type `P`.
+ * @param mr A memory resource used for allocating intermediate memory.
+ * @param vector The vector which to check for contiguity.
+ * @return true If the vector is contiguous on `P`.
+ * @return false Otherwise.
+ */
+template <TRACCC_CONSTRAINT(std::semiregular) P,
+          TRACCC_CONSTRAINT(std::equality_comparable) T>
+#if __cpp_concepts >= 201907L
+requires std::regular_invocable<P, T>
+#endif
+    bool is_contiguous_on(P&& projection, vecmem::data::vector_view<T> vector) {
+    // Grab the number of elements in our vector.
+    uint32_t n = vector.size();
+
+    // Get the output type of the projection.
+    using projection_t = std::invoke_result_t<P, T>;
+
+    // Allocate memory for intermediate values and outputs, then set them up.
+    std::unique_ptr<projection_t[]> iout = std::make_unique<projection_t[]>(n);
+    uint32_t iout_size = 0;
+
+    // Create a device vector
+    vecmem::device_vector<T> in(vector);
+
+    // Compress adjacent elements
+    for (std::size_t i = 0; i < n; ++i) {
+        if (i == 0) {
+            iout[iout_size++] = projection(in.at(i));
+        } else {
+            projection_t v = projection(in.at(i));
+
+            if (v != iout[iout_size - 1]) {
+                iout[iout_size++] = v;
+            }
+        }
+    }
+
+    // Check whether all elements are unique
+    std::unordered_set<projection_t> seen;
+
+    for (std::size_t i = 0; i < iout_size; ++i) {
+        projection_t& v = iout[i];
+
+        if (seen.count(v) == 1) {
+            return false;
+        } else {
+            seen.insert(v);
+        }
+    }
+
+    return true;
+}
+}  // namespace traccc::host

core/include/traccc/sanity/ordered_on.hpp

@@ -0,0 +1,68 @@
+/**
+ * traccc library, part of the ACTS project (R&D line)
+ *
+ * (c) 2024 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s).
+#include "traccc/definitions/concepts.hpp"
+
+// VecMem include(s).
+#include <vecmem/containers/data/vector_view.hpp>
+#include <vecmem/containers/device_vector.hpp>
+#include <vecmem/memory/memory_resource.hpp>
+#include <vecmem/utils/copy.hpp>
+
+// System include
+#if __cpp_concepts >= 201907L
+#include <concepts>
+#endif
+
+namespace traccc::host {
+/**
+ * @brief Sanity check that a given vector is ordered on a given relation.
+ *
+ * For a vector $v$ to be ordered on a relation $R$, it must be the case that
+ * for all indices $i$ and $j$, if $i < j$, then $R(i, j)$.
+ *
+ * @note This function runs in O(n) time.
+ *
+ * @note Although functions like `std::sort` requires the relation to be strict
+ * weak order, this function is more lax in its requirements. Rather, the
+ * relation should be a total preorder, i.e. a non-strict weak order.
+ *
+ * @note For any strict weak order $R$, `is_ordered_on(sort(R, v))` is true.
+ *
+ * @tparam R The type of relation $R$, a callable which returns a bool if the
+ * first argument can be immediately before the second type.
+ * @tparam T The type of the vector.
+ * @param relation A relation object of type `R`.
+ * @param mr A memory resource used for allocating intermediate memory.
+ * @param vector The vector which to check for ordering.
+ * @return true If the vector is ordered on `R`.
+ * @return false Otherwise.
+ */
+template <TRACCC_CONSTRAINT(std::semiregular) R, typename T>
+#if __cpp_concepts >= 201907L
+requires std::relation<R, T, T>
+#endif
+    bool is_ordered_on(R relation, vecmem::data::vector_view<T> vector) {
+    // Grab the number of elements in our vector.
+    uint32_t n = vector.size();
+
+    vecmem::device_vector<T> in(vector);
+
+    // Check for orderedness.
+    for (std::size_t i = 1; i < n; ++i) {
+        if (!relation(in.at(i - 1), in.at(i))) {
+            return false;
+        }
+    }
+
+    return true;
+}
+}  // namespace traccc::host

core/include/traccc/utils/projections.hpp

@@ -0,0 +1,29 @@
+/**
+ * traccc library, part of the ACTS project (R&D line)
+ *
+ * (c) 2024 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include "traccc/definitions/qualifiers.hpp"
+#include "traccc/edm/cell.hpp"
+#include "traccc/edm/measurement.hpp"
+
+namespace traccc {
+struct [[maybe_unused]] cell_module_projection{
+    TRACCC_HOST_DEVICE [[maybe_unused]] auto operator()(const traccc::cell& m)
+        const {return m.module_link;
+}
+}
+;
+
+struct [[maybe_unused]] measurement_module_projection{
+    TRACCC_HOST_DEVICE auto operator()(const traccc::measurement& m)
+        const {return m.module_link;
+}
+}
+;
+}  // namespace traccc

core/include/traccc/utils/relations.hpp

@@ -0,0 +1,34 @@
+/**
+ * traccc library, part of the ACTS project (R&D line)
+ *
+ * (c) 2024 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include "traccc/definitions/qualifiers.hpp"
+#include "traccc/edm/cell.hpp"
+
+namespace traccc {
+struct [[maybe_unused]] channel0_major_cell_order_relation{
+    TRACCC_HOST_DEVICE [[maybe_unused]] bool operator()(const traccc::cell& a,
+                                                        const traccc::cell& b)
+        const {if (a.module_link ==
+                   b.module_link){if (a.channel1 <= b.channel1){return true;
+}
+else if (a.channel1 == b.channel1) {
+    return a.channel0 <= b.channel0;
+}
+else {
+    return false;
+}
+}
+else {
+    return true;
+}
+}
+}
+;
+}  // namespace traccc

core/src/clusterization/sparse_ccl_algorithm.cpp

@@ -9,6 +9,10 @@
 #include "traccc/clusterization/sparse_ccl_algorithm.hpp"
 
 #include "traccc/clusterization/details/sparse_ccl.hpp"
+#include "traccc/sanity/contiguous_on.hpp"
+#include "traccc/sanity/ordered_on.hpp"
+#include "traccc/utils/projections.hpp"
+#include "traccc/utils/relations.hpp"
 
 // VecMem include(s).
 #include <vecmem/containers/device_vector.hpp>
@@ -22,6 +26,9 @@ sparse_ccl_algorithm::sparse_ccl_algorithm(vecmem::memory_resource& mr)
 sparse_ccl_algorithm::output_type sparse_ccl_algorithm::operator()(
     const cell_collection_types::const_view& cells_view) const {
 
+    assert(is_contiguous_on(cell_module_projection(), cells_view));
+    assert(is_ordered_on(channel0_major_cell_order_relation(), cells_view));
+
     // Run SparseCCL to fill CCL indices.
     const cell_collection_types::const_device cells{cells_view};
     vecmem::vector<unsigned int> cluster_indices{cells.size(), &(m_mr.get())};

device/cuda/CMakeLists.txt

@@ -28,6 +28,8 @@ traccc_add_library( traccc_cuda cuda TYPE SHARED
   "src/utils/opaque_stream.cpp"
   "src/utils/utils.hpp"
   "src/utils/utils.cpp"
+  "src/sanity/contiguous_on.cuh"
+  "src/sanity/ordered_on.cuh"
   # Seed finding code.
   "include/traccc/cuda/seeding/track_params_estimation.hpp"
   "src/seeding/track_params_estimation.cu"

device/cuda/src/clusterization/clusterization_algorithm.cu

@@ -6,10 +6,14 @@
  */
 
 // CUDA Library include(s).
+#include "../sanity/contiguous_on.cuh"
+#include "../sanity/ordered_on.cuh"
 #include "../utils/barrier.hpp"
 #include "../utils/cuda_error_handling.hpp"
 #include "../utils/utils.hpp"
 #include "traccc/cuda/clusterization/clusterization_algorithm.hpp"
+#include "traccc/utils/projections.hpp"
+#include "traccc/utils/relations.hpp"
 
 // Project include(s)
 #include "traccc/clusterization/device/ccl_kernel.hpp"
@@ -60,6 +64,11 @@ clusterization_algorithm::output_type clusterization_algorithm::operator()(
     const cell_collection_types::const_view& cells,
     const cell_module_collection_types::const_view& modules) const {
 
+    assert(is_contiguous_on(cell_module_projection(), m_mr.main, m_copy,
+                            m_stream, cells));
+    assert(is_ordered_on(channel0_major_cell_order_relation(), m_mr.main,
+                         m_copy, m_stream, cells));
+
     // Get a convenience variable for the stream that we'll be using.
     cudaStream_t stream = details::get_stream(m_stream);
 

device/cuda/src/finding/finding_algorithm.cu

@@ -6,10 +6,12 @@
  */
 
 // Project include(s).
+#include "../sanity/contiguous_on.cuh"
 #include "../utils/cuda_error_handling.hpp"
 #include "../utils/utils.hpp"
 #include "traccc/cuda/finding/finding_algorithm.hpp"
 #include "traccc/definitions/primitives.hpp"
+#include "traccc/definitions/qualifiers.hpp"
 #include "traccc/edm/device/finding_global_counter.hpp"
 #include "traccc/finding/candidate_link.hpp"
 #include "traccc/finding/device/add_links_for_holes.hpp"
@@ -20,6 +22,7 @@
 #include "traccc/finding/device/make_barcode_sequence.hpp"
 #include "traccc/finding/device/propagate_to_next_surface.hpp"
 #include "traccc/finding/device/prune_tracks.hpp"
+#include "traccc/utils/projections.hpp"
 
 // detray include(s).
 #include "detray/core/detector.hpp"
@@ -43,10 +46,10 @@
 #include <thrust/unique.h>
 
 // System include(s).
+#include <cassert>
 #include <vector>
 
 namespace traccc::cuda {
-
 namespace kernels {
 
 /// CUDA kernel for running @c traccc::device::make_barcode_sequence
@@ -270,6 +273,9 @@ finding_algorithm<stepper_t, navigator_t>::operator()(
     measurement_collection_types::const_view::size_type n_measurements =
         m_copy.get_size(measurements);
 
+    assert(is_contiguous_on(measurement_module_projection(), m_mr.main, m_copy,
+                            m_stream, measurements));
+
     // Get copy of barcode uniques
     measurement_collection_types::buffer uniques_buffer{n_measurements,
                                                         m_mr.main};