Merge pull request #964 from beomki-yeo/remove-full-jacobian

Remove full jacobian from stepper

Author: beomki-yeo

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

@@ -94,26 +94,12 @@ struct parameter_transporter : actor {
         // from curvilinear frame)
         if (!bound_params.surface_link().is_invalid()) {
 
-            // Previous surface
-            tracking_surface prev_sf{navigation.detector(),
-                                     bound_params.surface_link()};
-
-            const bound_to_free_matrix_t bound_to_free_jacobian =
-                prev_sf.bound_to_free_jacobian(gctx, bound_params);
-
-            auto vol = navigation.get_volume();
-            const auto vol_mat_ptr =
-                vol.has_material() ? vol.material_parameters(stepping().pos())
-                                   : nullptr;
-            stepping.set_full_jacobian(
-                sf.template visit_mask<get_full_jacobian_kernel>(
-                    sf.transform(gctx), bound_to_free_jacobian, vol_mat_ptr,
-                    propagation._stepping));
+            const auto full_jacobian = get_full_jacobian(propagation);
 
             // Calculate surface-to-surface covariance transport
-            const bound_matrix_t new_cov =
-                stepping.full_jacobian() * bound_params.covariance() *
-                matrix::transpose(stepping.full_jacobian());
+            const bound_matrix_t new_cov = full_jacobian *
+                                           bound_params.covariance() *
+                                           matrix::transpose(full_jacobian);
 
             stepping.bound_params().set_covariance(new_cov);
         }
@@ -126,6 +112,39 @@ struct parameter_transporter : actor {
         bound_params.set_surface_link(sf.barcode());
     }
 
+    template <typename propagator_state_t>
+    DETRAY_HOST_DEVICE inline bound_matrix_t get_full_jacobian(
+        propagator_state_t& propagation) const {
+
+        const auto& stepping = propagation._stepping;
+        const auto& navigation = propagation._navigation;
+
+        // Geometry context for this track
+        const auto& gctx = propagation._context;
+
+        // Current Surface
+        const auto sf = navigation.get_surface();
+
+        // Bound track params of departure surface
+        auto& bound_params = stepping.bound_params();
+
+        // Previous surface
+        tracking_surface prev_sf{navigation.detector(),
+                                 bound_params.surface_link()};
+
+        const bound_to_free_matrix_t bound_to_free_jacobian =
+            prev_sf.bound_to_free_jacobian(gctx, bound_params);
+
+        auto vol = navigation.get_volume();
+        const auto vol_mat_ptr = vol.has_material()
+                                     ? vol.material_parameters(stepping().pos())
+                                     : nullptr;
+
+        return sf.template visit_mask<get_full_jacobian_kernel>(
+            sf.transform(gctx), bound_to_free_jacobian, vol_mat_ptr,
+            propagation._stepping);
+    }
+
 };  // namespace detray
 
 }  // namespace detray

core/include/detray/propagator/base_stepper.hpp

@@ -181,18 +181,6 @@ class base_stepper {
             return m_jac_transport;
         }
 
-        /// @returns the current full Jacbian.
-        DETRAY_HOST_DEVICE
-        inline const bound_matrix_type &full_jacobian() const {
-            return m_full_jacobian;
-        }
-
-        /// Set new full Jacbian.
-        DETRAY_HOST_DEVICE
-        inline void set_full_jacobian(const bound_matrix_type &jac) {
-            m_full_jacobian = jac;
-        }
-
         /// Reset transport Jacbian.
         DETRAY_HOST_DEVICE
         inline void reset_transport_jacobian() {
@@ -230,10 +218,6 @@ class base_stepper {
         /// Jacobian transport matrix
         free_matrix_type m_jac_transport = matrix::identity<free_matrix_type>();
 
-        /// Full jacobian
-        bound_matrix_type m_full_jacobian =
-            matrix::identity<bound_matrix_type>();
-
         /// Bound covariance
         bound_track_parameters_type m_bound_params;
 

tests/integration_tests/cpu/propagator/jacobian_validation.cpp

@@ -438,7 +438,9 @@ struct bound_getter : actor {
             actor_state.m_path_length = stepping.path_length();
             actor_state.m_abs_path_length = stepping.abs_path_length();
             actor_state.m_param_destination = stepping.bound_params();
-            actor_state.m_jacobi = stepping.full_jacobian();
+            actor_state.m_jacobi =
+                parameter_transporter<algebra_t>().get_full_jacobian(
+                    propagation);
 
             // Stop navigation if the destination surface found
             propagation._heartbeat &= navigation.exit();