Add some extra assertions

Author: niermann999

core/include/detray/builders/detector_builder.hpp

@@ -69,6 +69,7 @@ class detector_builder {
     /// functionality
     template <class builder_t>
     DETRAY_HOST auto decorate(dindex volume_idx) -> builder_t* {
+        assert(has_volume(volume_idx));
 
         m_volumes[volume_idx] =
             std::make_unique<builder_t>(std::move(m_volumes[volume_idx]));

core/include/detray/builders/volume_builder.hpp

@@ -211,9 +211,8 @@ class volume_builder : public volume_builder_interface<detector_t> {
         std::size_t n_portals{0u};
         for (auto& sf_desc : m_surfaces) {
 
-            const auto sf = geometry::surface{det, sf_desc};
-
-            sf.template visit_mask<detail::mask_index_update>(sf_desc);
+            det._masks.template visit<detail::mask_index_update>(sf_desc.mask(),
+                                                                 sf_desc);
             sf_desc.set_volume(m_volume.index());
             sf_desc.update_transform(trf_offset);
             sf_desc.set_index(sf_offset++);
@@ -297,21 +296,6 @@ struct mask_index_update {
     }
 };
 
-/// TODO: Remove once the material builder is used everywhere
-/// A functor to update the material index in surface objects
-struct material_index_update {
-
-    template <typename group_t, typename index_t, typename surface_t>
-    DETRAY_HOST inline void operator()(
-        [[maybe_unused]] const group_t& group,
-        [[maybe_unused]] const index_t& /*index*/,
-        [[maybe_unused]] surface_t& sf) const {
-        if constexpr (!concepts::grid<typename group_t::value_type>) {
-            sf.update_material(static_cast<dindex>(group.size()));
-        }
-    }
-};
-
 }  // namespace detail
 
 }  // namespace detray

core/include/detray/geometry/surface.hpp

@@ -53,7 +53,11 @@ class surface {
     /// @param desc from that detector.
     DETRAY_HOST_DEVICE
     constexpr surface(const detector_t &det, const descr_t &desc)
-        : m_detector{det}, m_desc{desc} {}
+        : m_detector{det}, m_desc{desc} {
+        assert(!m_desc.barcode().is_invalid());
+        assert(m_desc.index() < det.surfaces().size());
+        assert(m_desc.transform() < det.transform_store().size());
+    }
 
     /// Constructor from detector @param det and barcode @param bcd in
     /// that detector.
@@ -84,20 +88,30 @@ class surface {
     /// @returns the surface barcode
     DETRAY_HOST_DEVICE
     constexpr auto barcode() const -> geometry::barcode {
+        assert(!m_desc.barcode().is_invalid());
         return m_desc.barcode();
     }
 
     /// @returns the index of the mother volume
     DETRAY_HOST_DEVICE
-    constexpr auto volume() const -> dindex { return barcode().volume(); }
+    constexpr auto volume() const -> dindex {
+        assert(barcode().volume() < m_detector.volumes().size());
+        return barcode().volume();
+    }
 
     /// @returns the index of the surface in the detector surface lookup
     DETRAY_HOST_DEVICE
-    constexpr auto index() const -> dindex { return barcode().index(); }
+    constexpr auto index() const -> dindex {
+        assert(barcode().index() < m_detector.surfaces().size());
+        return barcode().index();
+    }
 
     /// @returns the surface id (sensitive, passive or portal)
     DETRAY_HOST_DEVICE
-    constexpr auto id() const -> surface_id { return barcode().id(); }
+    constexpr auto id() const -> surface_id {
+        assert(barcode().id() != surface_id::e_unknown);
+        return barcode().id();
+    }
 
     /// @returns the extra bits in the barcode
     DETRAY_HOST_DEVICE
@@ -110,7 +124,7 @@ class surface {
     /// @returns the surface source link
     DETRAY_HOST_DEVICE
     constexpr auto source() const {
-        return m_detector.surface(m_desc.barcode()).source;
+        return m_detector.surface(barcode()).source;
     }
 
     /// @returns true if the surface is a senstive detector module.
@@ -156,6 +170,7 @@ class surface {
     DETRAY_HOST_DEVICE
     constexpr auto transform(const context &ctx) const
         -> const transform3_type & {
+        assert(m_desc.transform() < m_detector.transform_store().size());
         return m_detector.transform_store().at(m_desc.transform(), ctx);
     }
 
@@ -239,8 +254,8 @@ class surface {
     /// @tparam Args      types of additional arguments to the functor
     template <typename functor_t, typename... Args>
     DETRAY_HOST_DEVICE constexpr auto visit_mask(Args &&... args) const {
+        assert(!m_desc.mask().is_invalid());
         const auto &masks = m_detector.mask_store();
-
         return masks.template visit<functor_t>(m_desc.mask(),
                                                std::forward<Args>(args)...);
     }
@@ -251,8 +266,8 @@ class surface {
     /// @tparam Args      types of additional arguments to the functor
     template <typename functor_t, typename... Args>
     DETRAY_HOST_DEVICE constexpr auto visit_material(Args &&... args) const {
+        assert(has_material());
         const auto &materials = m_detector.material_store();
-
         return materials.template visit<functor_t>(m_desc.material(),
                                                    std::forward<Args>(args)...);
     }

core/include/detray/geometry/tracking_volume.hpp

@@ -57,7 +57,10 @@ class tracking_volume {
     /// Constructor from detector @param det and volume descriptor
     /// @param vol_idx from that detector.
     constexpr tracking_volume(const detector_t &det, const descr_t &desc)
-        : m_detector{det}, m_desc{desc} {}
+        : m_detector{det}, m_desc{desc} {
+        assert(m_desc.index() < det.volumes().size());
+        assert(m_desc.id() != volume_id::e_unknown);
+    }
 
     /// Constructor from detector @param det and volume index @param vol_idx in
     /// that detector.
@@ -168,8 +171,8 @@ class tracking_volume {
     /// @tparam Args      types of additional arguments to the functor
     template <typename functor_t, typename... Args>
     DETRAY_HOST_DEVICE constexpr auto visit_material(Args &&... args) const {
+        assert(has_material());
         const auto &materials = m_detector.material_store();
-
         return materials.template visit<functor_t>(m_desc.material(),
                                                    std::forward<Args>(args)...);
     }

core/include/detray/navigation/direct_navigator.hpp

@@ -288,6 +288,7 @@ class direct_navigator {
     DETRAY_HOST_DEVICE inline void init(const track_t &track, state &navigation,
                                         const navigation::config &cfg,
                                         const context_type &ctx) const {
+        assert(!track.is_invalid());
 
         if (navigation.is_complete()) {
             navigation.m_heartbeat = false;
@@ -305,6 +306,8 @@ class direct_navigator {
         const context_type &ctx = {},
         const bool is_before_actor_run = true) const {
 
+        assert(!track.is_invalid());
+
         if (navigation.is_complete()) {
             navigation.m_heartbeat = false;
             return true;

core/include/detray/navigation/navigator.hpp

@@ -521,6 +521,7 @@ class navigator {
         DETRAY_HOST_DEVICE
         inline void set_next(dindex pos) {
             m_next = pos;
+            assert(m_next >= 0);
             assert(m_next <= m_last + 1);
             assert(m_next < static_cast<dist_t>(k_cache_capacity) + 1);
         }
@@ -530,6 +531,7 @@ class navigator {
         inline void set_next(candidate_itr_t new_next) {
             m_next = static_cast<dist_t>(
                 detray::ranges::distance(m_candidates.begin(), new_next));
+            assert(m_next >= 0);
             assert(m_next < static_cast<dist_t>(k_cache_capacity));
         }
 
@@ -616,6 +618,7 @@ class navigator {
             const scalar_type mask_tol_scalor,
             const scalar_type overstep_tol) const {
 
+            assert(!sf_descr.barcode().is_invalid());
             const auto sf = geometry::surface{det, sf_descr};
 
             sf.template visit_mask<intersection_initialize<ray_intersector>>(
@@ -650,6 +653,8 @@ class navigator {
         const auto &det = navigation.detector();
         const auto volume = tracking_volume{det, navigation.volume()};
 
+        assert(!track.is_invalid());
+
         // Clean up state
         navigation.clear();
         navigation.m_heartbeat = true;
@@ -698,6 +703,9 @@ class navigator {
         const track_t &track, state &navigation, const navigation::config &cfg,
         const context_type &ctx = {},
         [[maybe_unused]] const bool is_before_actor = true) const {
+
+        assert(!track.is_invalid());
+
         // Candidates are re-evaluated based on the current trust level.
         // Should result in 'full trust'
         bool is_init = update_kernel(track, navigation, cfg, ctx);

core/include/detray/propagator/actors/parameter_resetter.hpp

@@ -35,6 +35,7 @@ struct parameter_resetter : actor {
         const auto sf = navigation.get_surface();
         stepping() = sf.bound_to_free_vector(propagation._context,
                                              stepping.bound_params());
+        assert(!stepping().is_invalid());
 
         // Reset jacobian transport to identity matrix
         stepping.reset_transport_jacobian();

core/include/detray/propagator/actors/parameter_transporter.hpp

@@ -111,6 +111,8 @@ struct parameter_transporter : actor {
 
         // Set surface link
         bound_params.set_surface_link(sf.barcode());
+
+        assert(!bound_params.is_invalid());
     }
 
     template <typename propagator_state_t>

core/include/detray/propagator/actors/pointwise_material_interactor.hpp

@@ -185,6 +185,8 @@ struct pointwise_material_interactor : actor {
                     interactor_state.projected_scattering_angle);
             }
         }
+
+        assert(!bound_params.is_invalid());
     }
 
     /// @brief Update the q over p of bound track parameter

core/include/detray/propagator/propagator.hpp

@@ -179,13 +179,15 @@ struct propagator {
         auto &stepping = propagation._stepping;
         auto &context = propagation._context;
         const auto &track = stepping();
+        assert(!track.is_invalid());
 
         // Initialize the navigation
         m_navigator.init(track, navigation, m_cfg.navigation, context);
         propagation._heartbeat = navigation.is_alive();
 
         // Run all registered actors/aborters after init
         run_actors(actor_state_refs, propagation);
+        assert(!track.is_invalid());
 
         // Find next candidate
         m_navigator.update(track, navigation, m_cfg.navigation, context);
@@ -210,6 +212,7 @@ struct propagator {
         auto &stepping = propagation._stepping;
         auto &context = propagation._context;
         const auto &track = stepping();
+        assert(!track.is_invalid());
 
         // Set access to the volume material for the stepper
         auto vol = navigation.get_volume();
@@ -234,6 +237,7 @@ struct propagator {
 
         // Run all registered actors/aborters after update
         run_actors(actor_state_refs, propagation);
+        assert(!track.is_invalid());
 
         // And check the status
         is_init |= m_navigator.update(track, navigation, m_cfg.navigation,
@@ -353,6 +357,7 @@ struct propagator {
                     break;
                 } else {
                     run_actors(actor_state_refs, propagation);
+                    assert(!track.is_invalid());
 
                     // And check the status
                     is_init |= m_navigator.update(
@@ -365,6 +370,7 @@ struct propagator {
 
                 // Synchronized actor
                 run_actors(actor_state_refs, propagation);
+                assert(!track.is_invalid());
 
                 // And check the status
                 is_init |= m_navigator.update(track, navigation,

core/include/detray/propagator/rk_stepper.ipp

@@ -415,6 +415,7 @@ DETRAY_HOST_DEVICE inline auto detray::rk_stepper<
         bvec2[1u] = bvec2_tmp[1u];
         bvec2[2u] = bvec2_tmp[2u];
 
+        assert(delta != 0.f);
         const vector3_type gradient = (bvec1 - bvec2) * (1.f / (2.f * delta));
 
         getter::element(dBdr, 0u, i) = gradient[0u];
@@ -475,6 +476,7 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper<
         return 0.f;
     }
 
+    assert(qop != 0.f);
     const scalar_type q = this->particle_hypothesis().charge();
     const scalar_type p = q / qop;
     const scalar_type mass = this->particle_hypothesis().mass();
@@ -487,6 +489,7 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper<
 
     // Assert that a momentum is a positive value
     assert(p >= 0.f);
+    assert(q != 0.f);
 
     // d(qop)ds, which is equal to (qop) * E * (-dE/ds) / p^2
     // or equal to (qop)^3 * E * (-dE/ds) / q^2
@@ -532,6 +535,8 @@ detray::rk_stepper<magnetic_field_t, algebra_t, constraint_t, policy_t,
 
     // Check Eq 3.12 of
     // (https://iopscience.iop.org/article/10.1088/1748-0221/4/04/P04016/meta)
+    assert(E2 != 0.f);
+    assert(g != 0.f);
     return dqopds * (1.f / qop * (3.f - p2 / E2) + 1.f / g * dgdqop);
 }
 
@@ -545,6 +550,10 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper<
          const detray::stepping::config& cfg, const bool do_reset,
          const material<scalar_type>* vol_mat_ptr) const {
 
+    // Check navigator and actor results
+    assert(dist_to_next != 0.f);
+    assert(!stepping().is_invalid());
+
     // Get stepper and navigator states
     auto& magnetic_field = stepping.m_magnetic_field;
 
@@ -564,6 +573,9 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper<
 
     // First Runge-Kutta point
     const auto bvec = magnetic_field.at(pos[0], pos[1], pos[2]);
+    assert(math::isfinite(bvec[0]));
+    assert(math::isfinite(bvec[1]));
+    assert(math::isfinite(bvec[2]));
     sd.b_first[0] = bvec[0];
     sd.b_first[1] = bvec[1];
     sd.b_first[2] = bvec[2];
@@ -577,6 +589,7 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper<
 
     /// RKN step trial and error estimation
     const auto estimate_error = [&](const scalar_type& h) {
+        assert(h != 0);
         // State the square and half of the step size
         const scalar_type h2{h * h};
         const scalar_type half_h{h * 0.5f};
@@ -587,6 +600,9 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper<
         const point3_type pos1 =
             pos + half_h * sd.t[0u] + h2 * 0.125f * sd.dtds[0u];
         const auto bvec1 = magnetic_field.at(pos1[0], pos1[1], pos1[2]);
+        assert(math::isfinite(bvec1[0]));
+        assert(math::isfinite(bvec1[1]));
+        assert(math::isfinite(bvec1[2]));
         sd.b_middle[0] = bvec1[0];
         sd.b_middle[1] = bvec1[1];
         sd.b_middle[2] = bvec1[2];
@@ -609,6 +625,9 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper<
         // Eq (84) of https://doi.org/10.1016/0029-554X(81)90063-X
         const point3_type pos2 = pos + h * sd.t[0u] + h2 * 0.5f * sd.dtds[2u];
         const auto bvec2 = magnetic_field.at(pos2[0], pos2[1], pos2[2]);
+        assert(math::isfinite(bvec2[0]));
+        assert(math::isfinite(bvec2[1]));
+        assert(math::isfinite(bvec2[2]));
         sd.b_last[0] = bvec2[0];
         sd.b_last[1] = bvec2[1];
         sd.b_last[2] = bvec2[2];
@@ -620,8 +639,7 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper<
 
         // Compute and check the local integration error estimate
         // @Todo
-        constexpr const scalar_type one_sixth{
-            static_cast<scalar_type>(1. / 6.)};
+        constexpr auto one_sixth{static_cast<scalar_type>(1. / 6.)};
         const vector3_type err_vec =
             one_sixth * h2 *
             (sd.dtds[0u] - sd.dtds[1u] - sd.dtds[2u] + sd.dtds[3u]);
@@ -633,6 +651,7 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper<
     /// estimate @param err
     const auto step_size_scaling =
         [&cfg](const scalar_type& err) -> scalar_type {
+        assert(err != 0.f);
         return static_cast<scalar_type>(
             math::min(math::max(math::sqrt(math::sqrt(cfg.rk_error_tol / err)),
                                 static_cast<scalar_type>(0.25)),
@@ -649,6 +668,7 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper<
 
         error = math::max(estimate_error(stepping.step_size()),
                           static_cast<scalar_type>(1e-20));
+        assert(math::isfinite(error));
 
         // Error is small enough
         // ---> break and advance track
@@ -683,6 +703,7 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper<
 
     // Advance track state
     stepping.advance_track(sd, vol_mat_ptr);
+    assert(!stepping().is_invalid());
 
     // Advance jacobian transport
     if (cfg.do_covariance_transport) {