diff --git a/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp b/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp index be8f2068e9..340ef10434 100644 --- a/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp +++ b/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp @@ -130,6 +130,9 @@ combinatorial_kalman_filter( std::vector> tips; + std::vector> sf_sequences; + sf_sequences.resize(config.max_track_candidates_per_track); + // Create propagator auto prop_cfg{config.propagation}; prop_cfg.navigation.estimate_scattering_noise = false; @@ -507,6 +510,10 @@ combinatorial_kalman_filter( typename detray::pathlimit_aborter::state aborter_state; + typename detray::surface_sequencer< + typename detector_t::surface_type>::state sequencer_state{ + vecmem::device_vector{ + vecmem::get_data(sf_sequences[link_id])}}; typename detray::parameter_transporter< typename detector_t::algebra_type>::state transporter_state; traccc::details::ckf_interactor_t::state interactor_state; @@ -533,7 +540,7 @@ combinatorial_kalman_filter( TRACCC_DEBUG_HOST("Propagating... "); propagator.propagate( propagation, - detray::tie(aborter_state, transporter_state, + detray::tie(aborter_state, sequencer_state, transporter_state, interaction_register_state, interactor_state, resetter_state, momentum_aborter_state, ckf_aborter_state)); diff --git a/core/include/traccc/finding/details/combinatorial_kalman_filter_types.hpp b/core/include/traccc/finding/details/combinatorial_kalman_filter_types.hpp index 089293b570..529870fe70 100644 --- a/core/include/traccc/finding/details/combinatorial_kalman_filter_types.hpp +++ b/core/include/traccc/finding/details/combinatorial_kalman_filter_types.hpp @@ -37,8 +37,10 @@ using ckf_interactor_t = detray::pointwise_material_interactor; /// Actor chain used in the Combinatorial Kalman Filter (CKF) +template using ckf_actor_chain_t = detray::actor_chain, + detray::surface_sequencer, detray::parameter_transporter, interaction_register, ckf_interactor_t, @@ -53,6 +55,6 @@ template using ckf_propagator_t = detray::propagator, detray::caching_navigator>, - ckf_actor_chain_t>; + ckf_actor_chain_t>; } // namespace traccc::details diff --git a/core/include/traccc/finding/finding_config.hpp b/core/include/traccc/finding/finding_config.hpp index 74b7a68fe7..5b82c8f5ce 100644 --- a/core/include/traccc/finding/finding_config.hpp +++ b/core/include/traccc/finding/finding_config.hpp @@ -10,6 +10,7 @@ // traccc include(s). #include "traccc/definitions/common.hpp" #include "traccc/definitions/primitives.hpp" +#include "traccc/fitting/fitting_config.hpp" #include "traccc/utils/particle.hpp" // detray include(s). @@ -46,6 +47,11 @@ struct finding_config { /// Enable the MBF smoother bool run_mbf_smoother = true; + /// Enable the Kalman smoother + bool run_kalman_smoother = false; + // Fitting config for the Kalman smoother + traccc::fitting_config kalman_smoother; + /// Minimum step length that track should make to reach the next surface. It /// should be set higher than the overstep tolerance not to make it stay on /// the same surface diff --git a/core/include/traccc/fitting/details/kalman_fitting.hpp b/core/include/traccc/fitting/details/kalman_fitting.hpp index ddd80dc3d0..970321fa69 100644 --- a/core/include/traccc/fitting/details/kalman_fitting.hpp +++ b/core/include/traccc/fitting/details/kalman_fitting.hpp @@ -71,7 +71,8 @@ typename edm::track_container::host kalman_fitting( seqs_buffer{ static_cast::size_type>( - std::max(fitted_track.constituent_links().size() * + std::max(static_cast( + fitted_track.constituent_links().size()) * fitter.config().surface_sequence_size_factor, fitter.config().min_surface_sequence_capacity)), mr, vecmem::data::buffer_type::resizable}; diff --git a/core/include/traccc/fitting/fitting_config.hpp b/core/include/traccc/fitting/fitting_config.hpp index 1aec19716f..5c01351f80 100644 --- a/core/include/traccc/fitting/fitting_config.hpp +++ b/core/include/traccc/fitting/fitting_config.hpp @@ -20,7 +20,7 @@ namespace traccc { /// Configuration struct for track fitting struct fitting_config { - std::size_t n_iterations = 1; + unsigned int n_iterations = 1u; /// Propagation configuration detray::propagation::config propagation{}; @@ -35,8 +35,8 @@ struct fitting_config { /// Smoothing with backward filter traccc::scalar covariance_inflation_factor = 1e3f; - std::size_t surface_sequence_size_factor = 5; - std::size_t min_surface_sequence_capacity = 100; + unsigned int surface_sequence_size_factor = 4u; + unsigned int min_surface_sequence_capacity = 100u; }; } // namespace traccc diff --git a/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp b/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp index affbb4f9ed..0aaf37608a 100644 --- a/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp +++ b/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp @@ -347,6 +347,7 @@ struct kalman_actor : detray::actor { auto& stepping = propagation._stepping; auto& navigation = propagation._navigation; + TRACCC_INFO_HOST_DEVICE("IN SMOOTHER"); TRACCC_VERBOSE_HOST_DEVICE("In Kalman actor (status %d)...", actor_state.fit_result); @@ -367,6 +368,7 @@ struct kalman_actor : detray::actor { // triggered only for sensitive surfaces if (navigation.is_on_sensitive()) { + TRACCC_INFO_HOST_DEVICE("On surface"); TRACCC_DEBUG_HOST( "-> on surface: " << navigation.current_surface()); @@ -379,9 +381,11 @@ struct kalman_actor : detray::actor { actor_state.add_to_sequence( std::as_const(navigation).current().sf_desc); } + DETRAY_ERROR_HOST_DEVICE("Surface not matched!"); return; } else if (actor_state.fit_result != kalman_fitter_status::SUCCESS) { + DETRAY_ERROR_HOST_DEVICE("ABORTED"); // Surface matched but encountered error: Abort fit navigation.abort(fitter_debug_msg{actor_state.fit_result}); propagation._heartbeat = false; @@ -435,7 +439,7 @@ struct kalman_actor : detray::actor { direction_e == kalman_actor_direction::BIDIRECTIONAL) { // Backward filter for smoothing - TRACCC_DEBUG_HOST_DEVICE("Run smoothing..."); + TRACCC_INFO_HOST_DEVICE("Run smoothing..."); // Forward filter did not find this state: cannot smoothe if (trk_state.filtered_params().is_invalid()) { @@ -449,6 +453,12 @@ struct kalman_actor : detray::actor { two_filters_smoother{}( trk_state, actor_state.m_measurements, bound_param, is_line); + if (!trk_state.is_smoothed()) { + TRACCC_ERROR_HOST_DEVICE( + "Track state NOT smoothed!"); + } else { + TRACCC_INFO_HOST_DEVICE("Track state IS smoothed!"); + } } } else { assert(false); @@ -487,6 +497,7 @@ struct kalman_actor : detray::actor { if (actor_state.finished() && !actor_state.do_precise_hole_count) { navigation.exit(); propagation._heartbeat = false; + TRACCC_INFO_HOST_DEVICE("Navigation exit"); return; } diff --git a/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp b/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp index 3c024e6fea..226a9f89d4 100644 --- a/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp +++ b/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp @@ -295,13 +295,14 @@ class kalman_fitter { /// @param fitter_state the state of kalman fitter [[nodiscard]] TRACCC_HOST_DEVICE kalman_fitter_status smooth(state& fitter_state) const { - TRACCC_VERBOSE_HOST_DEVICE("Run smoothing..."); + TRACCC_INFO_HOST_DEVICE("Run smoothing..."); if (fitter_state.m_fit_actor_state.sequencer().overflow()) { - TRACCC_ERROR_HOST_DEVICE("Barcode sequence overlow"); + TRACCC_ERROR_HOST_DEVICE("Surface sequence overlow"); return kalman_fitter_status::ERROR_BARCODE_SEQUENCE_OVERFLOW; } if (fitter_state.m_fit_actor_state.sequencer().sequence().empty()) { + TRACCC_ERROR_HOST_DEVICE("Surface sequence empty"); return kalman_fitter_status::ERROR_UPDATER_SKIPPED_STATE; } @@ -314,6 +315,12 @@ class kalman_fitter { while (!fitter_state.m_fit_actor_state.finished() && (!fitter_state.m_fit_actor_state.is_state() || fitter_state.m_fit_actor_state().is_hole())) { + if (!fitter_state.m_fit_actor_state.is_state()) { + TRACCC_ERROR_HOST_DEVICE("Not a track state"); + } + if (fitter_state.m_fit_actor_state().is_hole()) { + TRACCC_ERROR_HOST_DEVICE("HOLE"); + } fitter_state.m_fit_actor_state.next(); } @@ -367,6 +374,7 @@ class kalman_fitter { while (propagation._navigation.has_next_external() && propagation._navigation.next_external().barcode() != last.smoothed_params().surface_link()) { + TRACCC_INFO_HOST_DEVICE("Barcode does not match!"); TRACCC_DEBUG_HOST( "Advancing to next external surface from: " << propagation._navigation.next_external().barcode()); @@ -375,6 +383,7 @@ class kalman_fitter { // No valid states produced by forward pass if (propagation._navigation.finished()) { + TRACCC_ERROR_HOST_DEVICE("No matching barcodes!"); return kalman_fitter_status::ERROR_UPDATER_SKIPPED_STATE; } @@ -390,6 +399,7 @@ class kalman_fitter { m_cfg.covariance_inflation_factor); // Run the smoothing + TRACCC_INFO_HOST_DEVICE("Propagating"); propagator.propagate(propagation, fitter_state.backward_actor_state()); // Reset the backward mode to false @@ -398,10 +408,12 @@ class kalman_fitter { // Encountered error in the smoother? if (fitter_state.m_fit_actor_state.fit_result != kalman_fitter_status::SUCCESS) { + TRACCC_ERROR_HOST_DEVICE("Fit failed!"); return fitter_state.m_fit_actor_state.fit_result; } // Encountered error during propagation? if (!propagator.finished(propagation)) { + TRACCC_ERROR_HOST_DEVICE("-> Propagation failed!"); return kalman_fitter_status::ERROR_PROPAGATION_FAILURE; } @@ -519,6 +531,7 @@ class kalman_fitter { fit_res.constituent_links()) { assert(link.type == edm::track_constituent_link::track_state); auto trk_state = track_states.at(link.index); + // Fitting fails if any of non-hole track states is not smoothed if (!trk_state.is_hole() && !trk_state.is_smoothed()) { TRACCC_ERROR_HOST_DEVICE("Not all smoothed"); @@ -529,7 +542,6 @@ class kalman_fitter { } // Fitting succeeds if any of non-hole track states is not smoothed - TRACCC_DEBUG_HOST_DEVICE("Fit status: SUCCESS"); fit_res.fit_outcome() = track_fit_outcome::SUCCESS; return; } diff --git a/device/alpaka/src/finding/combinatorial_kalman_filter.hpp b/device/alpaka/src/finding/combinatorial_kalman_filter.hpp index e307ae999b..eaf4d68c60 100644 --- a/device/alpaka/src/finding/combinatorial_kalman_filter.hpp +++ b/device/alpaka/src/finding/combinatorial_kalman_filter.hpp @@ -150,7 +150,7 @@ struct propagate_to_next_surface { struct build_tracks { template ALPAKA_FN_ACC void operator()( - TAcc const& acc, bool run_mbf, + TAcc const& acc, bool run_mbf, const bool run_kalman_smoother, const device::build_tracks_payload payload) const { device::global_index_t globalThreadIdx = diff --git a/device/common/include/traccc/finding/device/build_tracks.hpp b/device/common/include/traccc/finding/device/build_tracks.hpp index 05c88f1086..23ca1831ae 100644 --- a/device/common/include/traccc/finding/device/build_tracks.hpp +++ b/device/common/include/traccc/finding/device/build_tracks.hpp @@ -64,7 +64,7 @@ struct build_tracks_payload { /// @param[inout] payload The function call payload /// TRACCC_HOST_DEVICE inline void build_tracks( - global_index_t globalIndex, bool run_mbf, + global_index_t globalIndex, bool run_mbf, const bool run_kalman_smoother, const build_tracks_payload& payload); } // namespace traccc::device diff --git a/device/common/include/traccc/finding/device/impl/build_tracks.ipp b/device/common/include/traccc/finding/device/impl/build_tracks.ipp index e24e274cda..d2e87b8275 100644 --- a/device/common/include/traccc/finding/device/impl/build_tracks.ipp +++ b/device/common/include/traccc/finding/device/impl/build_tracks.ipp @@ -20,7 +20,7 @@ namespace traccc::device { TRACCC_HOST_DEVICE inline void build_tracks( const global_index_t globalIndex, bool run_mbf, - const build_tracks_payload& payload) { + const bool run_kalman_smoother, const build_tracks_payload& payload) { const edm::measurement_collection::const_device measurements(payload.tracks_view.measurements); @@ -91,22 +91,33 @@ TRACCC_HOST_DEVICE inline void build_tracks( assert(L.meas_idx < n_meas); - if (run_mbf) { - const unsigned int track_state_index = + auto track_state_index{std::numeric_limits::max()}; + if (run_mbf || run_kalman_smoother) { + track_state_index = track_states.push_back(edm::make_track_state( measurements, L.meas_idx)); auto track_state = track_states.at(track_state_index); + // Get the filtered track parameters from the CKF + const auto& filtered_params = link_filtered_params.at(link_idx); + track_state.filtered_params() = filtered_params; + track_state.filtered_chi2() = L.chi2; track_state.set_hole(false); + } + + // Run the MBF smoother + if (run_mbf) { + auto track_state = track_states.at(track_state_index); + track_state.set_smoothed(true); // TODO: The fact that we store the chi2 three times is nonsense. - track_state.filtered_chi2() = L.chi2; - track_state.smoothed_chi2() = L.chi2; - track_state.backward_chi2() = L.chi2; + track_state.smoothed_chi2() = + std::numeric_limits::max(); + track_state.backward_chi2() = + std::numeric_limits::max(); const auto& predicted_params = link_predicted_params.at(link_idx); - const auto& filtered_params = link_filtered_params.at(link_idx); const auto& predicted_vec = predicted_params.vector(); const auto& predicted_covariance = predicted_params.covariance(); @@ -174,10 +185,10 @@ TRACCC_HOST_DEVICE inline void build_tracks( (predicted_covariance * big_lambda_tilde * predicted_covariance); - track_state.filtered_params() = filtered_params; - accumulated_jacobian = payload.jacobian_ptr[link_idx]; + *it = {edm::track_constituent_link::track_state, track_state_index}; + } else if (run_kalman_smoother) { *it = {edm::track_constituent_link::track_state, track_state_index}; } else { *it = {edm::track_constituent_link::measurement, L.meas_idx}; diff --git a/device/common/include/traccc/finding/device/impl/propagate_to_next_surface.ipp b/device/common/include/traccc/finding/device/impl/propagate_to_next_surface.ipp index 5e6a928b01..3b56ac5c28 100644 --- a/device/common/include/traccc/finding/device/impl/propagate_to_next_surface.ipp +++ b/device/common/include/traccc/finding/device/impl/propagate_to_next_surface.ipp @@ -84,21 +84,27 @@ TRACCC_HOST_DEVICE inline void propagate_to_next_surface( // Pathlimit aborter typename detray::detail::tuple_element<0, actor_tuple_type>::type::state s0{}; - typename detray::detail::tuple_element<1, actor_tuple_type>::type::state - s1{}; + // Surface sequencer + typename detray::detail::tuple_element<1, actor_tuple_type>::type::state s1{ + vecmem::device_vector< + typename propagator_t::detector_type::surface_type>( + *(payload.surfaces_view.ptr() + param_id))}; + // Parameter transporter state + typename detray::detail::tuple_element<2, actor_tuple_type>::type::state + s2{}; // CKF-interactor - typename detray::detail::tuple_element<3, actor_tuple_type>::type::state - s3{}; + typename detray::detail::tuple_element<4, actor_tuple_type>::type::state + s4{}; // Interaction register - typename detray::detail::tuple_element<2, actor_tuple_type>::type::state s2{ - s3}; + typename detray::detail::tuple_element<3, actor_tuple_type>::type::state s3{ + s4}; // Parameter resetter - typename detray::detail::tuple_element<4, actor_tuple_type>::type::state s4{ + typename detray::detail::tuple_element<5, actor_tuple_type>::type::state s5{ prop_cfg}; // Momentum aborter - typename detray::detail::tuple_element<5, actor_tuple_type>::type::state s5; - // CKF aborter typename detray::detail::tuple_element<6, actor_tuple_type>::type::state s6; + // CKF aborter + typename detray::detail::tuple_element<7, actor_tuple_type>::type::state s7; /* * If we are running the MBF smoother, we need to accumulate the Jacobians @@ -111,19 +117,23 @@ TRACCC_HOST_DEVICE inline void propagate_to_next_surface( payload.tmp_jacobian_ptr[param_id] = matrix::identity< bound_matrix>(); - s1._full_jacobian_ptr = &payload.tmp_jacobian_ptr[param_id]; + s2._full_jacobian_ptr = &payload.tmp_jacobian_ptr[param_id]; + } + + if (cfg.run_kalman_smoother) { } - s5.min_pT(static_cast(cfg.min_pT)); - s5.min_p(static_cast(cfg.min_p)); - s6.min_step_length = cfg.min_step_length_for_next_surface; - s6.max_count = cfg.max_step_counts_for_next_surface; + s6.min_pT(static_cast(cfg.min_pT)); + s6.min_p(static_cast(cfg.min_p)); + s7.min_step_length = cfg.min_step_length_for_next_surface; + s7.max_count = cfg.max_step_counts_for_next_surface; // Propagate to the next surface - propagator.propagate(propagation, detray::tie(s0, s1, s2, s3, s4, s5, s6)); + propagator.propagate(propagation, + detray::tie(s0, s1, s2, s3, s4, s5, s6, s7)); // If a surface found, add the parameter for the next step - if (s6.success) { + if (s7.success) { assert(propagation._navigation.is_on_sensitive()); assert(!propagation._stepping.bound_params().is_invalid()); diff --git a/device/common/include/traccc/finding/device/propagate_to_next_surface.hpp b/device/common/include/traccc/finding/device/propagate_to_next_surface.hpp index f9ca985819..d816a36595 100644 --- a/device/common/include/traccc/finding/device/propagate_to_next_surface.hpp +++ b/device/common/include/traccc/finding/device/propagate_to_next_surface.hpp @@ -82,6 +82,13 @@ struct propagate_to_next_surface_payload { bound_matrix* tmp_jacobian_ptr; + + /** + * @brief View object to the output barcode sequence + */ + vecmem::data::jagged_vector_view< + typename propagator_t::detector_type::surface_type> + surfaces_view; }; /// Function for propagating the kalman-updated tracks to the next surface diff --git a/device/common/include/traccc/fitting/device/fit_backward.hpp b/device/common/include/traccc/fitting/device/fit_backward.hpp index 7a11bdf0af..fe06614901 100644 --- a/device/common/include/traccc/fitting/device/fit_backward.hpp +++ b/device/common/include/traccc/fitting/device/fit_backward.hpp @@ -42,7 +42,7 @@ TRACCC_HOST_DEVICE inline void fit_backward( *(payload.surfaces_view.ptr() + param_id), fitter.config().propagation); - kalman_fitter_status fit_status = fitter.smooth(fitter_state); + const kalman_fitter_status fit_status = fitter.smooth(fitter_state); fitter.update_statistics(fitter_state); @@ -51,11 +51,11 @@ TRACCC_HOST_DEVICE inline void fit_backward( fitter.check_fitting_result(fitter_state, kalman_fitter_status::SUCCESS, fit_status); - if (fit_status == kalman_fitter_status::SUCCESS) { + /*if (fit_status == kalman_fitter_status::SUCCESS) { track = fitter_state.m_fit_res; } else { param_liveness.at(param_id) = 0u; - } + }*/ // TODO: Grab the smoothed state for next it } diff --git a/device/common/include/traccc/fitting/device/fit_prelude.hpp b/device/common/include/traccc/fitting/device/fit_prelude.hpp index 9f28acae4b..7a1e3e6601 100644 --- a/device/common/include/traccc/fitting/device/fit_prelude.hpp +++ b/device/common/include/traccc/fitting/device/fit_prelude.hpp @@ -35,15 +35,15 @@ TRACCC_HOST_DEVICE inline void fit_prelude( return; } - typename edm::track_state_collection::device track_states( - tracks_view.states); + // typename edm::track_state_collection::device track_states( + // tracks_view.states); vecmem::device_vector param_ids(param_ids_view); vecmem::device_vector param_liveness(param_liveness_view); const unsigned int param_id = param_ids.at(globalIndex); - auto track = tracks.at(param_id); + /*auto track = tracks.at(param_id); const typename edm::track_collection::const_device track_candidates{track_candidates_view.tracks}; @@ -62,7 +62,7 @@ TRACCC_HOST_DEVICE inline void fit_prelude( } // TODO: Set other stuff in the header? - track.params() = track_candidate.params(); + track.params() = track_candidate.params();*/ param_liveness.at(param_id) = 1u; } diff --git a/device/cuda/src/finding/combinatorial_kalman_filter.cuh b/device/cuda/src/finding/combinatorial_kalman_filter.cuh index 325d5b7702..45a359c311 100644 --- a/device/cuda/src/finding/combinatorial_kalman_filter.cuh +++ b/device/cuda/src/finding/combinatorial_kalman_filter.cuh @@ -8,6 +8,9 @@ #pragma once // Local include(s). +#include "../fitting/kernels/fill_fitting_sort_keys.hpp" +#include "../fitting/kernels/fit_backward.hpp" +#include "../fitting/kernels/fit_prelude.hpp" #include "../sanity/contiguous_on.cuh" #include "../utils/barrier.hpp" #include "../utils/cuda_error_handling.hpp" @@ -34,6 +37,8 @@ #include "traccc/finding/details/combinatorial_kalman_filter_types.hpp" #include "traccc/finding/device/barcode_surface_comparator.hpp" #include "traccc/finding/finding_config.hpp" +#include "traccc/fitting/details/kalman_fitting_types.hpp" +#include "traccc/fitting/device/fill_fitting_sort_keys.hpp" #include "traccc/utils/logging.hpp" #include "traccc/utils/memory_resource.hpp" #include "traccc/utils/projections.hpp" @@ -170,6 +175,10 @@ combinatorial_kalman_filter( bound_track_parameters_collection_types::buffer link_filtered_parameter_buffer(0, mr.main); + vecmem::data::jagged_vector_buffer + seqs_buffer{std::vector{0u}, mr.main, mr.host, + vecmem::data::buffer_type::resizable}; + /* * If we are aiming to run the MBF smoother at the end of the track * finding, we need some space to store the intermediate Jacobians @@ -182,11 +191,29 @@ combinatorial_kalman_filter( link_predicted_parameter_buffer = bound_track_parameters_collection_types::buffer( link_buffer_capacity, mr.main); + } + + if (config.run_mbf_smoother || config.run_kalman_smoother) { link_filtered_parameter_buffer = bound_track_parameters_collection_types::buffer( link_buffer_capacity, mr.main); } + // if (config.run_kalman_smoother) { + const unsigned int n_prop_streams{ + math::min(config.max_num_branches_per_seed * n_seeds, 200000u)}; + const unsigned int n_surfaces_per_track{ + std::max(config.max_track_candidates_per_track * + config.kalman_smoother.surface_sequence_size_factor, + config.kalman_smoother.min_surface_sequence_capacity)}; + std::vector seqs_sizes(n_prop_streams, n_surfaces_per_track); + + seqs_buffer = + vecmem::data::jagged_vector_buffer{ + seqs_sizes, mr.main, mr.host, vecmem::data::buffer_type::resizable}; + copy.setup(seqs_buffer)->ignore(); + //} + // Create a buffer of tip links vecmem::data::vector_buffer tips_buffer{ config.max_num_branches_per_seed * n_seeds, mr.main, @@ -271,6 +298,21 @@ combinatorial_kalman_filter( links_buffer = std::move(new_links_buffer); + if (config.run_mbf_smoother || config.run_kalman_smoother) { + bound_track_parameters_collection_types::buffer + new_link_filtered_parameter_buffer{link_buffer_capacity, + mr.main}; + + copy(link_filtered_parameter_buffer, + new_link_filtered_parameter_buffer) + ->wait(); + + TRACCC_CUDA_ERROR_CHECK(cudaStreamSynchronize(stream)); + + link_filtered_parameter_buffer = + std::move(new_link_filtered_parameter_buffer); + } + if (config.run_mbf_smoother) { vecmem::unique_alloc_ptr< bound_matrix[]> @@ -280,9 +322,6 @@ combinatorial_kalman_filter( bound_track_parameters_collection_types::buffer new_link_predicted_parameter_buffer{link_buffer_capacity, mr.main}; - bound_track_parameters_collection_types::buffer - new_link_filtered_parameter_buffer{link_buffer_capacity, - mr.main}; TRACCC_CUDA_ERROR_CHECK(cudaMemcpyAsync( new_jacobian_ptr.get(), jacobian_ptr.get(), @@ -293,17 +332,12 @@ combinatorial_kalman_filter( copy(link_predicted_parameter_buffer, new_link_predicted_parameter_buffer) ->wait(); - copy(link_filtered_parameter_buffer, - new_link_filtered_parameter_buffer) - ->wait(); TRACCC_CUDA_ERROR_CHECK(cudaStreamSynchronize(stream)); jacobian_ptr = std::move(new_jacobian_ptr); link_predicted_parameter_buffer = std::move(new_link_predicted_parameter_buffer); - link_filtered_parameter_buffer = - std::move(new_link_filtered_parameter_buffer); } } @@ -501,7 +535,8 @@ combinatorial_kalman_filter( .n_in_params = n_candidates, .tips_view = tips_buffer, .tip_lengths_view = tip_length_buffer, - .tmp_jacobian_ptr = tmp_jacobian_ptr.get()}; + .tmp_jacobian_ptr = tmp_jacobian_ptr.get(), + .surfaces_view = seqs_buffer}; const unsigned int nThreads = warp_size * 4; const unsigned int nBlocks = @@ -649,7 +684,7 @@ combinatorial_kalman_filter( unsigned int n_states; - if (config.run_mbf_smoother) { + if (config.run_mbf_smoother || config.run_kalman_smoother) { n_states = std::accumulate(tips_length_host.begin(), tips_length_host.end(), 0u); } else { @@ -682,12 +717,95 @@ combinatorial_kalman_filter( .link_filtered_parameter_view = link_filtered_parameter_buffer, }; kernels::build_tracks<<>>( - config.run_mbf_smoother, payload); + config.run_mbf_smoother, config.run_kalman_smoother, payload); TRACCC_CUDA_ERROR_CHECK(cudaGetLastError()); str.synchronize(); } + /// Run a Kalman filter in the backpropagation to smoothe the tracks + if (config.run_kalman_smoother) { + + std::cout << "KALMAN SMOOTHER" << std::endl; + + const typename edm::track_container:: + const_view track_candidates_view{track_candidates_buffer}; + + // Get the number of tracks. + const unsigned int n_tracks = + copy.get_size(track_candidates_view.tracks); + + // Create the result buffer. + /*typename edm::track_container::buffer track_states_buffer{ + {tips_length_host, mr.main, mr.host, + vecmem::data::buffer_type::resizable}, + {n_states, mr.main, vecmem::data::buffer_type::resizable}, + track_candidates_view.measurements}; + copy.setup(track_states_buffer.tracks)->ignore(); + copy.setup(track_states_buffer.states)->ignore(); + + // Return early, if there are no tracks. + if (n_tracks == 0) { + return track_states_buffer; + }*/ + + // Create the buffers for sorting the parameter IDs. + /*vecmem::data::vector_buffer keys_buffer(n_tracks, + mr.main); + vecmem::data::vector_buffer param_ids_buffer(n_tracks, + mr.main); + vecmem::data::vector_buffer param_liveness_buffer( + n_tracks, mr.main); + vecmem::copy::event_type keys_setup_event = copy.setup(keys_buffer); + vecmem::copy::event_type param_ids_setup_event = + copy.setup(param_ids_buffer); + vecmem::copy::event_type param_liveness_setup_event = + copy.setup(param_liveness_buffer); + keys_setup_event->ignore(); + param_ids_setup_event->ignore(); + param_liveness_setup_event->ignore(); + + // Launch parameters for all the kernels. + const unsigned int nThreads = warp_size * 4; + const unsigned int nBlocks = (n_tracks + nThreads - 1) / nThreads; + + // Fill the keys and param_ids buffers. + fill_fitting_sort_keys(nBlocks, nThreads, stream, + track_candidates_view.tracks, keys_buffer, + param_ids_buffer); + + // Sort the key to get the sorted parameter ids + vecmem::device_vector keys_device(keys_buffer); + vecmem::device_vector param_ids_device(param_ids_buffer); + thrust::sort_by_key( + thrust::cuda::par_nosync(std::pmr::polymorphic_allocator(&mr.main)) + .on(stream), + keys_device.begin(), keys_device.end(), param_ids_device.begin()); + + // Run the fitting, using the sorted parameter IDs. + fit_prelude(nBlocks, nThreads, 0, stream, param_ids_buffer, + track_candidates_view, track_candidates_buffer, + param_liveness_buffer); + TRACCC_CUDA_ERROR_CHECK(cudaGetLastError()); + str.synchronize();*/ + + // Allocate the fitting kernels's payload in host memory. + using fitter_t = traccc::details::kalman_fitter_t; + device::fit_payload host_payload{ + .det_data = det, + .field_data = field, + .param_ids_view = param_ids_buffer, + .param_liveness_view = param_liveness_buffer, + .tracks_view = track_candidates_buffer, + .surfaces_view = seqs_buffer}; + + // Run the track fitting + fit_backward(nBlocks, nThreads, 0, stream, + config.kalman_smoother, host_payload); + TRACCC_CUDA_ERROR_CHECK(cudaGetLastError()); + } + return track_candidates_buffer; } diff --git a/device/cuda/src/finding/kernels/build_tracks.cu b/device/cuda/src/finding/kernels/build_tracks.cu index 5d2607ca6c..473e35fa63 100644 --- a/device/cuda/src/finding/kernels/build_tracks.cu +++ b/device/cuda/src/finding/kernels/build_tracks.cu @@ -17,9 +17,10 @@ namespace traccc::cuda::kernels { -__global__ void build_tracks(bool run_mbf, +__global__ void build_tracks(bool run_mbf, const bool run_kalman_smoother, device::build_tracks_payload payload) { - device::build_tracks(details::global_index1(), run_mbf, payload); + device::build_tracks(details::global_index1(), run_mbf, run_kalman_smoother, + payload); } } // namespace traccc::cuda::kernels diff --git a/device/cuda/src/finding/kernels/build_tracks.cuh b/device/cuda/src/finding/kernels/build_tracks.cuh index cc21675b14..f603b13f86 100644 --- a/device/cuda/src/finding/kernels/build_tracks.cuh +++ b/device/cuda/src/finding/kernels/build_tracks.cuh @@ -13,6 +13,6 @@ namespace traccc::cuda::kernels { -__global__ void build_tracks(bool run_mbf, +__global__ void build_tracks(bool run_mbf, const bool run_kalman_smoother, device::build_tracks_payload payload); } diff --git a/device/sycl/src/finding/combinatorial_kalman_filter.hpp b/device/sycl/src/finding/combinatorial_kalman_filter.hpp index 3c5f57dab0..6b7e8f7da1 100644 --- a/device/sycl/src/finding/combinatorial_kalman_filter.hpp +++ b/device/sycl/src/finding/combinatorial_kalman_filter.hpp @@ -583,6 +583,7 @@ combinatorial_kalman_filter( ::sycl::nd_item<1> item) { device::build_tracks(details::global_index(item), false && config.run_mbf_smoother, + config.run_kalman_smoother, {.seeds_view = seeds, .links_view = links, .tips_view = tips, diff --git a/examples/options/src/track_finding.cpp b/examples/options/src/track_finding.cpp index efaea66eda..34450c9b21 100644 --- a/examples/options/src/track_finding.cpp +++ b/examples/options/src/track_finding.cpp @@ -93,6 +93,14 @@ track_finding::track_finding() : interface("Track Finding Options") { po::value(&m_config.duplicate_removal_minimum_length) ->default_value(m_config.duplicate_removal_minimum_length), "Minimum track length for deduplication (0 to disable) [cardinal]"); + m_desc.add_options()("finding-run-mbf-smoother", + po::value(&m_config.run_mbf_smoother) + ->default_value(m_config.run_mbf_smoother), + "Enable the MBF smoother"); + m_desc.add_options()("finding-run-kalman-smoother", + po::value(&m_config.run_kalman_smoother) + ->default_value(m_config.run_kalman_smoother), + "Enable the Kalman smoother"); } void track_finding::read(const po::variables_map &) { @@ -152,6 +160,11 @@ std::unique_ptr track_finding::as_printable() const { cat->add_child(std::make_unique( "Minimum p", std::format("{} GeV", m_config.min_p / traccc::unit::GeV))); + cat->add_child(std::make_unique( + "Run MBF smoother", m_config.run_mbf_smoother ? "true" : "false")); + cat->add_child(std::make_unique( + "Run Kalman smoother", + m_config.run_kalman_smoother ? "true" : "false")); return cat; }