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

ref: Update to detray version 0.85.0 and adapt to new algebra-plugins API #798

Merged
merged 1 commit into from
Dec 8, 2024
Merged

ref: Update to detray version 0.85.0 and adapt to new algebra-plugins API #798

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
8 changes: 5 additions & 3 deletions core/include/traccc/definitions/primitives.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
#elif ALGEBRA_PLUGINS_INCLUDE_SMATRIX
#include "traccc/plugins/algebra/smatrix_definitions.hpp"
#elif ALGEBRA_PLUGINS_INCLUDE_VC
#include "traccc/plugins/algebra/vc_definitions.hpp"
#include "traccc/plugins/algebra/vc_aos_definitions.hpp"
#elif ALGEBRA_PLUGINS_INCLUDE_VECMEM
#include "traccc/plugins/algebra/vecmem_definitions.hpp"
#endif
Expand All @@ -40,7 +40,9 @@ using point3 = detray::dpoint3D<default_algebra>;
using vector3 = detray::dvector3D<default_algebra>;
using variance3 = point3;
using transform3 = detray::dtransform3D<default_algebra>;
template <std::size_t ROW, std::size_t COL>
using matrix = detray::dmatrix<default_algebra, ROW, COL>;

namespace getter = detray::getter;
namespace vector = detray::vector;
namespace matrix = detray::matrix;

} // namespace traccc
37 changes: 15 additions & 22 deletions core/include/traccc/edm/track_state.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#pragma once

// Project include(s).
#include "traccc/definitions/primitives.hpp"
#include "traccc/definitions/qualifiers.hpp"
#include "traccc/edm/container.hpp"
#include "traccc/edm/measurement.hpp"
Expand Down Expand Up @@ -43,12 +44,11 @@ template <typename algebra_t>
struct track_state {

using scalar_type = detray::dscalar<algebra_t>;
using size_type = detray::dsize_type<algebra_t>;

using bound_track_parameters_type =
detray::bound_track_parameters<algebra_t>;
using bound_matrix = detray::bound_matrix<algebra_t>;
using matrix_operator = detray::dmatrix_operator<algebra_t>;
using size_type = detray::dsize_type<algebra_t>;
template <size_type ROWS, size_type COLS>
using matrix_type = detray::dmatrix<algebra_t, ROWS, COLS>;

Expand Down Expand Up @@ -85,17 +85,15 @@ struct track_state {
matrix_type<D, 1> ret;
switch (m_measurement.subs.get_indices()[0]) {
case e_bound_loc0:
matrix_operator().element(ret, 0, 0) = m_measurement.local[0];
getter::element(ret, 0, 0) = m_measurement.local[0];
if constexpr (D == 2u) {
matrix_operator().element(ret, 1, 0) =
m_measurement.local[1];
getter::element(ret, 1, 0) = m_measurement.local[1];
}
break;
case e_bound_loc1:
matrix_operator().element(ret, 0, 0) = m_measurement.local[1];
getter::element(ret, 0, 0) = m_measurement.local[1];
if constexpr (D == 2u) {
matrix_operator().element(ret, 1, 0) =
m_measurement.local[0];
getter::element(ret, 1, 0) = m_measurement.local[0];
}
break;
default:
Expand All @@ -115,23 +113,19 @@ struct track_state {
matrix_type<D, D> ret;
switch (m_measurement.subs.get_indices()[0]) {
case e_bound_loc0:
matrix_operator().element(ret, 0, 0) =
m_measurement.variance[0];
getter::element(ret, 0, 0) = m_measurement.variance[0];
if constexpr (D == 2u) {
matrix_operator().element(ret, 0, 1) = 0.f;
matrix_operator().element(ret, 1, 0) = 0.f;
matrix_operator().element(ret, 1, 1) =
m_measurement.variance[1];
getter::element(ret, 0, 1) = 0.f;
getter::element(ret, 1, 0) = 0.f;
getter::element(ret, 1, 1) = m_measurement.variance[1];
}
break;
case e_bound_loc1:
matrix_operator().element(ret, 0, 0) =
m_measurement.variance[1];
getter::element(ret, 0, 0) = m_measurement.variance[1];
if constexpr (D == 2u) {
matrix_operator().element(ret, 0, 1) = 0.f;
matrix_operator().element(ret, 1, 0) = 0.f;
matrix_operator().element(ret, 1, 1) =
m_measurement.variance[0];
getter::element(ret, 0, 1) = 0.f;
getter::element(ret, 1, 0) = 0.f;
getter::element(ret, 1, 1) = m_measurement.variance[0];
}
break;
default:
Expand Down Expand Up @@ -200,8 +194,7 @@ struct track_state {
private:
detray::geometry::barcode m_surface_link;
measurement m_measurement;
bound_matrix m_jacobian =
matrix_operator().template zero<e_bound_size, e_bound_size>();
bound_matrix m_jacobian = matrix::zero<bound_matrix>();
bound_track_parameters_type m_predicted;
scalar_type m_filtered_chi2 = 0.f;
bound_track_parameters_type m_filtered;
Expand Down
58 changes: 26 additions & 32 deletions core/include/traccc/fitting/kalman_filter/gain_matrix_smoother.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,11 @@ struct gain_matrix_smoother {

// Type declarations
using scalar_type = detray::dscalar<algebra_t>;
using matrix_operator = detray::dmatrix_operator<algebra_t>;
using size_type = detray::dsize_type<algebra_t>;
template <size_type ROWS, size_type COLS>
using matrix_type = detray::dmatrix<algebra_t, ROWS, COLS>;
using bound_vector_type = detray::bound_vector<algebra_t>;
using bound_matrix_type = detray::bound_matrix<algebra_t>;

/// Gain matrix smoother operation
///
Expand Down Expand Up @@ -69,44 +70,38 @@ struct gain_matrix_smoother {
const auto& cur_filtered = cur_state.filtered();

// Next track state parameters
const matrix_type<e_bound_size, e_bound_size>& next_jacobian =
next_state.jacobian();
const matrix_type<e_bound_size, 1>& next_smoothed_vec =
next_smoothed.vector();
const matrix_type<e_bound_size, e_bound_size>& next_smoothed_cov =
next_smoothed.covariance();
const matrix_type<e_bound_size, 1>& next_predicted_vec =
next_predicted.vector();
const matrix_type<e_bound_size, e_bound_size>& next_predicted_cov =
const bound_matrix_type& next_jacobian = next_state.jacobian();
const bound_vector_type& next_smoothed_vec = next_smoothed.vector();
const bound_matrix_type& next_smoothed_cov = next_smoothed.covariance();
const bound_vector_type& next_predicted_vec = next_predicted.vector();
const bound_matrix_type& next_predicted_cov =
next_predicted.covariance();

// Current track state parameters
const matrix_type<e_bound_size, 1>& cur_filtered_vec =
cur_filtered.vector();
const matrix_type<e_bound_size, e_bound_size>& cur_filtered_cov =
cur_filtered.covariance();
const bound_vector_type& cur_filtered_vec = cur_filtered.vector();
const bound_matrix_type& cur_filtered_cov = cur_filtered.covariance();

// Regularization matrix for numerical stability
static constexpr scalar_type epsilon = 1e-13f;
const matrix_type<e_bound_size, e_bound_size> regularization =
matrix_operator().template identity<e_bound_size, e_bound_size>() *
epsilon;
const matrix_type<e_bound_size, e_bound_size>
regularized_predicted_cov = next_predicted_cov + regularization;
constexpr scalar_type epsilon = 1e-13f;
const auto regularization =
matrix::identity<bound_matrix_type>() * epsilon;
const bound_matrix_type regularized_predicted_cov =
next_predicted_cov + regularization;

// Calculate smoothed parameter for current state
const matrix_type<e_bound_size, e_bound_size> A =
cur_filtered_cov * matrix_operator().transpose(next_jacobian) *
matrix_operator().inverse(regularized_predicted_cov);
const bound_matrix_type A = cur_filtered_cov *
matrix::transpose(next_jacobian) *
matrix::inverse(regularized_predicted_cov);

const matrix_type<e_bound_size, 1> smt_vec =
const bound_vector_type smt_vec =
cur_filtered_vec + A * (next_smoothed_vec - next_predicted_vec);
const matrix_type<e_bound_size, e_bound_size> smt_cov =
cur_filtered_cov + A * (next_smoothed_cov - next_predicted_cov) *
matrix_operator().transpose(A);
const bound_matrix_type smt_cov =
cur_filtered_cov +
A * (next_smoothed_cov - next_predicted_cov) * matrix::transpose(A);

cur_state.smoothed().set_vector(smt_vec);
cur_state.smoothed().set_covariance(smt_cov);

// Wrap the phi in the range of [-pi, pi]
wrap_phi(cur_state.smoothed());

Expand All @@ -127,12 +122,11 @@ struct gain_matrix_smoother {
const matrix_type<D, D>& V =
cur_state.template measurement_covariance<D>();
const matrix_type<D, 1> residual = meas_local - H * smt_vec;
const matrix_type<D, D> R =
V - H * smt_cov * matrix_operator().transpose(H);
const matrix_type<1, 1> chi2 = matrix_operator().transpose(residual) *
matrix_operator().inverse(R) * residual;
const matrix_type<D, D> R = V - H * smt_cov * matrix::transpose(H);
const matrix_type<1, 1> chi2 =
matrix::transpose(residual) * matrix::inverse(R) * residual;

cur_state.smoothed_chi2() = matrix_operator().element(chi2, 0, 0);
cur_state.smoothed_chi2() = getter::element(chi2, 0, 0);

return;
}
Expand Down
31 changes: 13 additions & 18 deletions core/include/traccc/fitting/kalman_filter/gain_matrix_updater.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,10 +19,11 @@ template <typename algebra_t>
struct gain_matrix_updater {

// Type declarations
using matrix_operator = detray::dmatrix_operator<algebra_t>;
using size_type = detray::dsize_type<algebra_t>;
template <size_type ROWS, size_type COLS>
using matrix_type = detray::dmatrix<algebra_t, ROWS, COLS>;
using bound_vector_type = detray::bound_vector<algebra_t>;
using bound_matrix_type = detray::bound_matrix<algebra_t>;

/// Gain matrix updater operation
///
Expand Down Expand Up @@ -70,12 +71,9 @@ struct gain_matrix_updater {
const auto meas = trk_state.get_measurement();

// Some identity matrices
// @Note: Make constexpr work
const matrix_type<e_bound_size, e_bound_size> I66 =
matrix_operator().template identity<e_bound_size, e_bound_size>();

const matrix_type<D, D> I_m =
matrix_operator().template identity<D, D>();
// @TODO: Make constexpr work
const auto I66 = matrix::identity<bound_matrix_type>();
const auto I_m = matrix::identity<matrix_type<D, D>>();

matrix_type<D, e_bound_size> H = meas.subs.template projector<D>();

Expand All @@ -84,12 +82,10 @@ struct gain_matrix_updater {
trk_state.template measurement_local<D>();

// Predicted vector of bound track parameters
const matrix_type<e_bound_size, 1>& predicted_vec =
bound_params.vector();
const bound_vector_type& predicted_vec = bound_params.vector();

// Predicted covaraince of bound track parameters
const matrix_type<e_bound_size, e_bound_size>& predicted_cov =
bound_params.covariance();
const bound_matrix_type& predicted_cov = bound_params.covariance();

// Set track state parameters
trk_state.predicted().set_vector(predicted_vec);
Expand All @@ -108,12 +104,11 @@ struct gain_matrix_updater {
trk_state.template measurement_covariance<D>();

const matrix_type<D, D> M =
H * predicted_cov * matrix_operator().transpose(H) + V;
H * predicted_cov * matrix::transpose(H) + V;

// Kalman gain matrix
const matrix_type<6, D> K = predicted_cov *
matrix_operator().transpose(H) *
matrix_operator().inverse(M);
const matrix_type<6, D> K =
predicted_cov * matrix::transpose(H) * matrix::inverse(M);

// Calculate the filtered track parameters
const matrix_type<6, 1> filtered_vec =
Expand All @@ -125,8 +120,8 @@ struct gain_matrix_updater {

// Calculate the chi square
const matrix_type<D, D> R = (I_m - H * K) * V;
const matrix_type<1, 1> chi2 = matrix_operator().transpose(residual) *
matrix_operator().inverse(R) * residual;
const matrix_type<1, 1> chi2 =
matrix::transpose(residual) * matrix::inverse(R) * residual;

// Return false if track is parallel to z-axis or phi is not finite
const scalar theta = bound_params.theta();
Expand All @@ -138,7 +133,7 @@ struct gain_matrix_updater {
// Set the track state parameters
trk_state.filtered().set_vector(filtered_vec);
trk_state.filtered().set_covariance(filtered_cov);
trk_state.filtered_chi2() = matrix_operator().element(chi2, 0, 0);
trk_state.filtered_chi2() = getter::element(chi2, 0, 0);

// Wrap the phi in the range of [-pi, pi]
wrap_phi(trk_state.filtered());
Expand Down
2 changes: 1 addition & 1 deletion core/include/traccc/seeding/detail/seeding_config.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,7 @@ struct seedfinder_config {

TRACCC_HOST_DEVICE
size_t get_num_rbins() const {
return static_cast<size_t>(rMax + getter::norm(beamPos));
return static_cast<size_t>(rMax + vector::norm(beamPos));
}

TRACCC_HOST_DEVICE
Expand Down
2 changes: 1 addition & 1 deletion core/include/traccc/seeding/spacepoint_binning_helper.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ inline TRACCC_HOST_DEVICE size_t is_valid_sp(const seedfinder_config& config,
if (spPhi > config.phiMax || spPhi < config.phiMin) {
return detray::detail::invalid_value<size_t>();
}
size_t r_index = static_cast<size_t>(getter::perp(
size_t r_index = static_cast<size_t>(vector::perp(
vector2{sp.x() - config.beamPos[0], sp.y() - config.beamPos[1]}));

if (r_index < config.get_num_rbins()) {
Expand Down
14 changes: 7 additions & 7 deletions core/include/traccc/seeding/track_params_estimation_helper.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,22 +86,22 @@ seed_to_bound_vector(const spacepoint_collection_t& sp_collection,
scalar B = uv2[1] - A * uv2[0];

// Radius (with a sign)
scalar R = -getter::perp(vector2{1.f, A}) / (2.f * B);
scalar R = -vector::perp(vector2{1.f, A}) / (2.f * B);
// The (1/tanTheta) of momentum in the new frame
scalar invTanTheta =
local2[2] / (2.f * R * math::asin(getter::perp(local2) / (2.f * R)));
local2[2] / (2.f * R * math::asin(vector::perp(local2) / (2.f * R)));

// The momentum direction in the new frame (the center of the circle
// has the coordinate (-1.*A/(2*B), 1./(2*B)))
vector3 transDirection =
vector3({1.f, A, scalar(getter::perp(vector2{1.f, A})) * invTanTheta});
vector3({1.f, A, scalar(vector::perp(vector2{1.f, A})) * invTanTheta});
// Transform it back to the original frame
vector3 direction =
transform3::rotate(trans._data, vector::normalize(transDirection));

// The estimated phi and theta
getter::element(params, e_bound_phi, 0) = getter::phi(direction);
getter::element(params, e_bound_theta, 0) = getter::theta(direction);
getter::element(params, e_bound_phi, 0) = vector::phi(direction);
getter::element(params, e_bound_theta, 0) = vector::theta(direction);

// The measured loc0 and loc1
const auto& meas_for_spB = spB.meas;
Expand All @@ -110,10 +110,10 @@ seed_to_bound_vector(const spacepoint_collection_t& sp_collection,

// The estimated q/pt in [GeV/c]^-1 (note that the pt is the
// projection of momentum on the transverse plane of the new frame)
scalar qOverPt = 1.f / (R * getter::norm(bfield));
scalar qOverPt = 1.f / (R * vector::norm(bfield));
// The estimated q/p in [GeV/c]^-1
getter::element(params, e_bound_qoverp, 0) =
qOverPt / getter::perp(vector2{1.f, invTanTheta});
qOverPt / vector::perp(vector2{1.f, invTanTheta});

// Make sure the time is a finite value
assert(std::isfinite(getter::element(params, e_bound_time, 0)));
Expand Down
7 changes: 2 additions & 5 deletions core/include/traccc/utils/seed_generator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,6 @@ namespace traccc {
template <typename detector_t>
struct seed_generator {
using algebra_type = typename detector_t::algebra_type;
using matrix_operator = detray::dmatrix_operator<algebra_type>;
using cxt_t = typename detector_t::geometry_context;

/// Constructor with detector
Expand Down Expand Up @@ -58,9 +57,7 @@ struct seed_generator {

const cxt_t ctx{};
auto bound_vec = sf.free_to_bound_vector(ctx, free_param);

auto bound_cov =
matrix_operator().template zero<e_bound_size, e_bound_size>();
auto bound_cov = matrix::zero<detray::bound_matrix<algebra_type>>();

bound_track_parameters bound_param{surface_link, bound_vec, bound_cov};

Expand All @@ -82,7 +79,7 @@ struct seed_generator {
bound_param[i] = std::normal_distribution<scalar>(
bound_param[i], m_stddevs[i])(m_generator);

matrix_operator().element(bound_param.covariance(), i, i) =
getter::element(bound_param.covariance(), i, i) =
m_stddevs[i] * m_stddevs[i];
}

Expand Down
Loading
Loading