Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ref: Remove more state data from the steppers and streamline the parameter resetter #879

Draft
wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

ref: Remove more state data from the steppers and streamline the parameter resetter #879

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions core/include/detray/navigation/policies.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ struct stepper_default_policy : actor {

// Not a severe change to track state expected
// Policy is called after stepsize update -> use prev. step size
if (math::fabs(stepping.prev_step_size()) <
if (math::fabs(stepping.step_size()) <
math::fabs(
stepping.constraints().template size<>(stepping.direction())) -
pol_state.tol) {
Expand Down Expand Up @@ -112,7 +112,7 @@ struct stepper_rk_policy : actor {
auto &navigation = propagation._navigation;

// Policy is called after stepsize update -> use prev. step size
const scalar rel_correction{(stepping.prev_step_size() - navigation()) /
const scalar rel_correction{(stepping.step_size() - navigation()) /
navigation()};

// Large correction to the stepsize - re-initialize the volume
Expand Down
57 changes: 49 additions & 8 deletions core/include/detray/propagator/actor_chain.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,19 +1,24 @@
/** Detray library, part of the ACTS project (R&D line)
*
* (c) 2022 CERN for the benefit of the ACTS project
* (c) 2022-2024 CERN for the benefit of the ACTS project
*
* Mozilla Public License Version 2.0
*/

#pragma once

#include <type_traits>
#include <utility>

// Project include(s)
#include "detray/definitions/detail/containers.hpp"
#include "detray/definitions/detail/qualifiers.hpp"
#include "detray/propagator/base_actor.hpp"
#include "detray/utils/tuple_helpers.hpp"

// System include(s)
#include <concepts>
#include <optional>
#include <type_traits>
#include <utility>

namespace detray {

/// The interface to the actors and aborters in the propagation.
Expand Down Expand Up @@ -44,7 +49,29 @@ class actor_chain {
}

/// @returns the actor list
DETRAY_HOST_DEVICE const actor_list_type &actors() const { return _actors; }
DETRAY_HOST_DEVICE const actor_list_type &actors() const {
return m_actors;
}

/// @returns a tuple of default constructible actor states and a
/// corresponding tuple of references
DETRAY_HOST_DEVICE
static constexpr auto make_actor_states() {
// Only possible if each state is default initializable
if constexpr ((std::default_initializable<typename actors_t::state> &&
...)) {
return tuple_t<typename actors_t::state...>{};
} else {
return std::nullopt;
}
}

/// @returns a tuple of reference for every state in the tuple @param t
DETRAY_HOST_DEVICE static constexpr state make_ref_tuple(
tuple_t<typename actors_t::state...> &t) {
return make_ref_tuple(t,
std::make_index_sequence<sizeof...(actors_t)>{});
}

private:
/// Call the actors. Either single actor or composition.
Expand All @@ -58,7 +85,13 @@ class actor_chain {
actor_states_t &states,
propagator_state_t &p_state) const {
if constexpr (!typename actor_t::is_comp_actor()) {
actr(detail::get<typename actor_t::state &>(states), p_state);
// No actor state defined (empty)
if constexpr (std::same_as<typename actor_t::state,
detray::actor::state>) {
actr(p_state);
} else {
actr(detail::get<typename actor_t::state &>(states), p_state);
}
} else {
actr(states, p_state);
}
Expand All @@ -73,11 +106,19 @@ class actor_chain {
DETRAY_HOST_DEVICE inline void run(
actor_states_t &states, propagator_state_t &p_state,
std::index_sequence<indices...> /*ids*/) const {
(run(detail::get<indices>(_actors), states, p_state), ...);
(run(detail::get<indices>(m_actors), states, p_state), ...);
}

/// @returns a tuple of reference for every state in the tuple @param t
template <std::size_t... indices>
DETRAY_HOST_DEVICE static constexpr state make_ref_tuple(
tuple_t<typename actors_t::state...> &t,
std::index_sequence<indices...> /*ids*/) {
return detray::tie(detail::get<indices>(t)...);
}

/// Tuple of actors
actor_list_type _actors = {};
actor_list_type m_actors = {};
};

/// Empty actor chain (placeholder)
Expand Down
24 changes: 0 additions & 24 deletions core/include/detray/propagator/actors/aborters.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,28 +95,4 @@ struct target_aborter : actor {
}
};

/// Aborter triggered when the next surface is reached
struct next_surface_aborter : actor {
struct state {
// minimal step length to prevent from staying on the same surface
scalar min_step_length = 0.f;
bool success = false;
};

template <typename propagator_state_t>
DETRAY_HOST_DEVICE void operator()(state &abrt_state,
propagator_state_t &prop_state) const {

auto &navigation = prop_state._navigation;
auto &stepping = prop_state._stepping;

// Abort at the next sensitive surface
if (navigation.is_on_sensitive() &&
stepping.path_from_surface() > abrt_state.min_step_length) {
prop_state._heartbeat &= navigation.abort();
abrt_state.success = true;
}
}
};

} // namespace detray
51 changes: 7 additions & 44 deletions core/include/detray/propagator/actors/parameter_resetter.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
/** Detray library, part of the ACTS project (R&D line)
*
* (c) 2022-2023 CERN for the benefit of the ACTS project
* (c) 2022-2024 CERN for the benefit of the ACTS project
*
* Mozilla Public License Version 2.0
*/
Expand All @@ -19,45 +19,8 @@ namespace detray {
template <typename algebra_t>
struct parameter_resetter : actor {

using scalar_type = dscalar<algebra_t>;

struct state {};

/// Mask store visitor
struct kernel {

// Matrix actor
using transform3_type = dtransform3D<algebra_t>;
using matrix_operator = dmatrix_operator<algebra_t>;

template <typename mask_group_t, typename index_t,
typename stepper_state_t>
DETRAY_HOST_DEVICE inline void operator()(
const mask_group_t& mask_group, const index_t& index,
const transform3_type& trf3, const dindex sf_idx,
stepper_state_t& stepping) const {

// Note: How is it possible with "range"???
const auto& mask = mask_group[index];

// Reset the free vector
stepping() = detail::bound_to_free_vector(trf3, mask,
stepping.bound_params());

// Reset the path length
stepping.reset_path_from_surface();

// Reset jacobian transport to identity matrix
stepping.reset_transport_jacobian();

// Reset the surface index
stepping.set_prev_sf_index(sf_idx);
}
};

template <typename propagator_state_t>
DETRAY_HOST_DEVICE void operator()(state& /*resetter_state*/,
propagator_state_t& propagation) const {
DETRAY_HOST_DEVICE void operator()(propagator_state_t& propagation) const {

const auto& navigation = propagation._navigation;
auto& stepping = propagation._stepping;
Expand All @@ -68,14 +31,14 @@ struct parameter_resetter : actor {
return;
}

using geo_cxt_t =
typename propagator_state_t::detector_type::geometry_context;
const geo_cxt_t ctx{};
typename propagator_state_t::detector_type::geometry_context ctx{};

// Surface
// Update free params after bound params were changed by actors
const auto sf = navigation.get_surface();
stepping() = sf.bound_to_free_vector(ctx, stepping.bound_params());

sf.template visit_mask<kernel>(sf.transform(ctx), sf.index(), stepping);
// Reset jacobian transport to identity matrix
stepping.reset_transport_jacobian();
}
};

Expand Down
77 changes: 34 additions & 43 deletions core/include/detray/propagator/actors/parameter_transporter.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,31 +21,27 @@ struct parameter_transporter : actor {

/// @name Type definitions for the struct
/// @{

using scalar_type = dscalar<algebra_t>;
// Transformation matching this struct
using transform3_type = dtransform3D<algebra_t>;
// scalar_type
using scalar_type = dscalar<algebra_t>;
// Matrix actor
using matrix_operator = dmatrix_operator<algebra_t>;
// 2D matrix type
template <std::size_t ROWS, std::size_t COLS>
using matrix_type = dmatrix<algebra_t, ROWS, COLS>;
// bound matrix type
using bound_matrix_t = bound_matrix<algebra_t>;
// Matrix type for bound to free jacobian
using bound_to_free_matrix_t = bound_to_free_matrix<algebra_t>;
/// @}

struct state {};

struct get_full_jacobian_kernel {

template <typename mask_group_t, typename index_t,
typename propagator_state_t>
typename stepper_state_t>
DETRAY_HOST_DEVICE inline bound_matrix_t operator()(
const mask_group_t& /*mask_group*/, const index_t& /*index*/,
const transform3_type& trf3,
const bound_to_free_matrix<algebra_t>& bound_to_free_jacobian,
const propagator_state_t& propagation) const {
const bound_to_free_matrix_t& bound_to_free_jacobian,
const material<scalar_type>* vol_mat_ptr,
const stepper_state_t& stepping) const {

using frame_t = typename mask_group_t::value_type::shape::
template local_frame_type<algebra_t>;
Expand All @@ -56,38 +52,28 @@ struct parameter_transporter : actor {
using free_to_bound_matrix_t =
typename jacobian_engine_t::free_to_bound_matrix_type;

// Stepper and Navigator states
auto& stepping = propagation._stepping;

// Free vector
const auto& free_params = stepping();

// Free to bound jacobian at the destination surface
const free_to_bound_matrix_t free_to_bound_jacobian =
jacobian_engine_t::free_to_bound_jacobian(trf3, free_params);

// Transport jacobian in free coordinate
const free_matrix_t& free_transport_jacobian =
stepping.transport_jacobian();
jacobian_engine_t::free_to_bound_jacobian(trf3, stepping());

// Path correction factor
free_matrix_t path_correction = jacobian_engine_t::path_correction(
stepping().pos(), stepping().dir(), stepping.dtds(),
stepping.dqopds(), trf3);
const free_matrix_t path_correction =
jacobian_engine_t::path_correction(
stepping().pos(), stepping().dir(), stepping.dtds(),
stepping.dqopds(vol_mat_ptr), trf3);

const free_matrix_t correction_term =
matrix_operator()
.template identity<e_free_size, e_free_size>() +
path_correction;

return free_to_bound_jacobian * correction_term *
free_transport_jacobian * bound_to_free_jacobian;
stepping.transport_jacobian() * bound_to_free_jacobian;
}
};

template <typename propagator_state_t>
DETRAY_HOST_DEVICE void operator()(state& /*actor_state*/,
propagator_state_t& propagation) const {
DETRAY_HOST_DEVICE void operator()(propagator_state_t& propagation) const {
auto& stepping = propagation._stepping;
const auto& navigation = propagation._navigation;

Expand All @@ -97,46 +83,51 @@ struct parameter_transporter : actor {
return;
}

using detector_type = typename propagator_state_t::detector_type;
using geo_cxt_t = typename detector_type::geometry_context;
const geo_cxt_t ctx{};
typename propagator_state_t::detector_type::geometry_context ctx{};

// Current Surface
const auto sf = navigation.get_surface();

// Bound track params of departure surface
auto& bound_params = stepping.bound_params();

// Covariance is transported only when the previous surface is an
// actual tracking surface. (i.e. This disables the covariance transport
// from curvilinear frame)
if (!detail::is_invalid_value(stepping.prev_sf_index())) {
if (!bound_params.surface_link().is_invalid()) {

// Previous surface
tracking_surface<detector_type> prev_sf{navigation.detector(),
stepping.prev_sf_index()};
tracking_surface prev_sf{navigation.detector(),
bound_params.surface_link()};

const bound_to_free_matrix<algebra_t> bound_to_free_jacobian =
prev_sf.bound_to_free_jacobian(ctx, stepping.bound_params());
const bound_to_free_matrix_t bound_to_free_jacobian =
prev_sf.bound_to_free_jacobian(ctx, 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(ctx), bound_to_free_jacobian, propagation));
sf.transform(ctx), bound_to_free_jacobian, vol_mat_ptr,
propagation._stepping));

// Calculate surface-to-surface covariance transport
const bound_matrix_t new_cov =
stepping.full_jacobian() *
stepping.bound_params().covariance() *
stepping.full_jacobian() * bound_params.covariance() *
matrix_operator().transpose(stepping.full_jacobian());

stepping.bound_params().set_covariance(new_cov);
}

// Convert free to bound vector
stepping.bound_params().set_parameter_vector(
bound_params.set_parameter_vector(
sf.free_to_bound_vector(ctx, stepping()));

// Set surface link
stepping.bound_params().set_surface_link(sf.barcode());

return;
bound_params.set_surface_link(sf.barcode());
}

}; // namespace detray

} // namespace detray
1 change: 0 additions & 1 deletion core/include/detray/propagator/actors/pointwise_material_interactor.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,6 @@ struct pointwise_material_interactor : actor {

struct state {

/// @TODO: Consider using the particle information in stepping::config
/// Evaluated energy loss
scalar_type e_loss{0.f};
/// Evaluated projected scattering angle
Expand Down
Loading
Loading