Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-05-27 07:24:20 

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008 
0009 // Project include(s)
0010 // clang-format off
0011 #include "detray/definitions/algebra.hpp"
0012 #include "detray/algebra/utils/approximately_equal.hpp"
0013 #include "detray/algebra/utils/casts.hpp"
0014 #include "detray/algebra/utils/print.hpp"
0015 // clang-format on
0016 
0017 // Test include(s)
0018 #include "detray/test/framework/types.hpp"
0019 
0020 // GoogleTest include(s).
0021 #include <gtest/gtest.h>
0022 
0023 // System include(s)
0024 #include <array>
0025 #include <limits>
0026 
0027 using namespace detray;
0028 
0029 using value_t = float;
0030 
0031 constexpr value_t tol{1e-5f};
0032 
0033 /// This test the vector functions on an SoA (Vc::Vector) based vector
0034 TEST(detray_algebra_soa, vector) {
0035   using test_algebra_t = detray::vc_soa<value_t>;
0036   using vector3_v = dvector3D<test_algebra_t>;
0037 
0038   // Value type is Vc::Vector<float>
0039   using scalar_t = dscalar<test_algebra_t>;
0040 
0041   static_assert(detray::concepts::scalar<scalar_t>);
0042   static_assert(detray::concepts::vector<vector3_v>);
0043   static_assert(detray::concepts::vector3D<vector3_v>);
0044 
0045   // Cast simd scalar to different precisions
0046   using scalar_f = Vc::float_v;
0047   using scalar_d = Vc::double_v;
0048   using scalar_i = Vc::int_v;
0049 
0050   auto s1 = detray::algebra::cast_to<float>(scalar_t(1.f));
0051   auto s2 = detray::algebra::cast_to<double>(scalar_t(2.f));
0052   auto s3 = detray::algebra::cast_to<int>(scalar_t(3.f));
0053 
0054   static_assert(std::same_as<decltype(s1), scalar_f>);
0055   static_assert(std::same_as<decltype(s2), scalar_d>);
0056   static_assert(std::same_as<decltype(s3), scalar_i>);
0057 
0058   ASSERT_TRUE((s1 == scalar_f(1.f)).isFull());
0059   ASSERT_TRUE((s2 == scalar_d(2.f)).isFull());
0060   ASSERT_TRUE((s3 == scalar_i(3.f)).isFull());
0061 
0062   vector3_v a{1.f, 2.f, 3.f};
0063   vector3_v b{4.f, 5.f, 6.f};
0064 
0065   EXPECT_TRUE((a[0] == scalar_t(1.f)).isFull());
0066   EXPECT_TRUE((a[1] == scalar_t(2.f)).isFull());
0067   EXPECT_TRUE((a[2] == scalar_t(3.f)).isFull());
0068 
0069   // Test printing
0070   std::cout << a << std::endl;
0071 
0072   // Test comparison
0073   constexpr auto epsilon{std::numeric_limits<value_t>::epsilon()};
0074 
0075   EXPECT_TRUE(detray::algebra::approx_equal(a, a));
0076   EXPECT_TRUE(detray::algebra::approx_equal(a, a, epsilon));
0077   EXPECT_FALSE(detray::algebra::approx_equal(a, b));
0078 
0079   value_t rel_err = 1.f + 10.f * epsilon;
0080   vector3_v a_err = rel_err * a;
0081   EXPECT_TRUE(detray::algebra::approx_equal(a, a_err, 11.f * epsilon));
0082   EXPECT_FALSE(detray::algebra::approx_equal(a, a_err, 9.f * epsilon));
0083 
0084   rel_err = 1.f + 17.f * epsilon;
0085   a_err = rel_err * a;
0086   EXPECT_TRUE(detray::algebra::approx_equal(a, a_err, 18.f * epsilon));
0087   EXPECT_FALSE(detray::algebra::approx_equal(a, a_err, 16.f * epsilon));
0088 
0089   // Swap an element
0090   vector3_v a_err_cpy = a_err;
0091   EXPECT_TRUE(a_err_cpy == a_err);
0092   EXPECT_TRUE(detray::algebra::approx_equal(a_err_cpy, a_err));
0093 
0094   auto& vec_elem = a_err[0];
0095   vec_elem[0] += 1.f;
0096   EXPECT_FALSE(a_err_cpy == a_err);
0097   EXPECT_FALSE(detray::algebra::approx_equal(a_err_cpy, a_err));
0098   // Cast simd vectors to different precision
0099   auto a_cast_f = detray::algebra::cast_to<float>(a);
0100   auto a_cast_d = detray::algebra::cast_to<double>(a);
0101   auto a_cast_i = detray::algebra::cast_to<int>(a);
0102 
0103   using algebra_f_t = detray::vc_soa<float>;
0104   using algebra_d_t = detray::vc_soa<double>;
0105   using algebra_i_t = detray::vc_soa<int>;
0106   static_assert(std::same_as<decltype(a_cast_f), dvector3D<algebra_f_t>>);
0107   static_assert(std::same_as<decltype(a_cast_d), dvector3D<algebra_d_t>>);
0108   static_assert(std::same_as<decltype(a_cast_i), dvector3D<algebra_i_t>>);
0109 
0110   for (int i = 0; i < 3; ++i) {
0111     EXPECT_TRUE(
0112         (a_cast_f[i] == detray::algebra::cast_to<float>(a[i])).isFull());
0113     EXPECT_TRUE(
0114         (a_cast_d[i] == detray::algebra::cast_to<double>(a[i])).isFull());
0115     EXPECT_TRUE((a_cast_i[i] == detray::algebra::cast_to<int>(a[i])).isFull());
0116   }
0117 
0118   // Masked comparison
0119   auto m = a.compare(a);
0120   EXPECT_TRUE(m[0].isFull());
0121   EXPECT_TRUE(m[1].isFull());
0122   EXPECT_TRUE(m[2].isFull());
0123 
0124   m = a.compare(b);
0125   EXPECT_FALSE(m[0].isFull());
0126   EXPECT_FALSE(m[1].isFull());
0127   EXPECT_FALSE(m[2].isFull());
0128 
0129   // Full comparisons
0130   EXPECT_TRUE(a == a);
0131   EXPECT_FALSE(a == b);
0132 
0133   // Addition
0134   auto v_add = a + b;
0135   EXPECT_TRUE((v_add[0] == scalar_t(5.f)).isFull());
0136   EXPECT_TRUE((v_add[1] == scalar_t(7.f)).isFull());
0137   EXPECT_TRUE((v_add[2] == scalar_t(9.f)).isFull());
0138 
0139   // Subration
0140   auto v_sub = a - b;
0141   EXPECT_TRUE((v_sub[0] == scalar_t(-3.f)).isFull());
0142   EXPECT_TRUE((v_sub[1] == scalar_t(-3.f)).isFull());
0143   EXPECT_TRUE((v_sub[2] == scalar_t(-3.f)).isFull());
0144 
0145   // Multiplication
0146   auto v_mul = a * b;
0147   EXPECT_TRUE((v_mul[0] == scalar_t(4.f)).isFull());
0148   EXPECT_TRUE((v_mul[1] == scalar_t(10.f)).isFull());
0149   EXPECT_TRUE((v_mul[2] == scalar_t(18.f)).isFull());
0150 
0151   // Division
0152   auto v_div = a / b;
0153   EXPECT_TRUE((v_div[0] == scalar_t(0.25f)).isFull());
0154   EXPECT_TRUE((v_div[1] == scalar_t(0.4f)).isFull());
0155   EXPECT_TRUE((v_div[2] == scalar_t(0.5f)).isFull());
0156 
0157   // Scalar multiplication
0158   auto v_smul = 2.f * b;
0159   EXPECT_TRUE((v_smul[0] == scalar_t(8.f)).isFull());
0160   EXPECT_TRUE((v_smul[1] == scalar_t(10.f)).isFull());
0161   EXPECT_TRUE((v_smul[2] == scalar_t(12.f)).isFull());
0162 
0163   // Expression
0164   auto v_expr = (b / a) - (2.5f * b) + vector3_v{};
0165   EXPECT_TRUE((v_expr[0] == scalar_t(-6.f)).isFull());
0166   EXPECT_TRUE((v_expr[1] == scalar_t(-10.f)).isFull());
0167   EXPECT_TRUE((v_expr[2] == scalar_t(-13.f)).isFull());
0168 
0169   auto d{vector::dot(a, b)};
0170   EXPECT_TRUE((d == scalar_t(32.f)).isFull());
0171 
0172   scalar_t norms_a{vector::norm(vector::normalize(a))};
0173   scalar_t norms_b{vector::norm(vector::normalize(b))};
0174   for (unsigned int i{0u}; i < norms_a.size(); ++i) {
0175     EXPECT_NEAR(norms_a[i], 1.f, tol);
0176     EXPECT_NEAR(norms_b[i], 1.f, tol);
0177   }
0178 
0179   auto cr{vector::cross(a, b)};
0180   EXPECT_TRUE((cr[0] == scalar_t(-3.f)).isFull());
0181   EXPECT_TRUE((cr[1] == scalar_t(6.f)).isFull());
0182   EXPECT_TRUE((cr[2] == scalar_t(-3.f)).isFull());
0183 
0184   static_assert(std::is_convertible_v<decltype(v_expr), vector3_v>,
0185                 "expression type not convertible");
0186 }
0187 
0188 /// This test the getter functions on an SoA (Vc::Vector) based vector
0189 TEST(detray_algebra_soa, getter) {
0190   using test_algebra_t = detray::vc_soa<value_t>;
0191 
0192   using vector3_v = dvector3D<test_algebra_t>;
0193 
0194   vector3_v a{1.f, 2.f, 3.f};
0195 
0196   // All results in the vector are the same, so only check the first one
0197 
0198   // Phi angle
0199   auto v_phi = vector::phi(a);
0200   EXPECT_NEAR(v_phi[0], static_cast<value_t>(std::atan2(2., 1.)), tol);
0201 
0202   // Perpendicular projection
0203   auto v_perp = vector::perp(a);
0204   EXPECT_NEAR(v_perp[0], std::sqrt(5.), tol);
0205 
0206   // Theta angle
0207   auto v_theta = vector::theta(a);
0208   EXPECT_NEAR(v_theta[0], static_cast<value_t>(std::atan2(std::sqrt(5.), 3.)),
0209               tol);
0210 
0211   // Norm of the vector
0212   auto v_norm = vector::norm(a);
0213   EXPECT_NEAR(v_norm[0], std::sqrt(14.), tol);
0214 
0215   // Eta of the vector
0216   auto v_eta = vector::eta(a);
0217   EXPECT_NEAR(v_eta[0],
0218               static_cast<value_t>(std::atanh(1. / std::sqrt(14.) * 3.)), tol);
0219 }
0220 
0221 /// This test an SoA (Vc::Vector) based affine transform3
0222 TEST(detray_algebra_soa, transform3) {
0223   using test_algebra_t = detray::vc_soa<value_t>;
0224 
0225   // Print the linear algebra types of this backend
0226   using algebra::operator<<;
0227 
0228   using vector3 = dvector3D<test_algebra_t>;
0229   using point3 = dpoint3D<test_algebra_t>;
0230   using scalar_t = dscalar<test_algebra_t>;
0231   using transform3 = dtransform3D<test_algebra_t>;
0232 
0233   static_assert(detray::concepts::transform3D<transform3>);
0234 
0235   transform3 idty{};
0236 
0237   EXPECT_TRUE((idty(0, 0) == scalar_t::One()).isFull());
0238   EXPECT_TRUE((idty(1, 0) == scalar_t::Zero()).isFull());
0239   EXPECT_TRUE((idty(2, 0) == scalar_t::Zero()).isFull());
0240   EXPECT_TRUE((idty(0, 1) == scalar_t::Zero()).isFull());
0241   EXPECT_TRUE((idty(1, 1) == scalar_t::One()).isFull());
0242   EXPECT_TRUE((idty(2, 1) == scalar_t::Zero()).isFull());
0243   EXPECT_TRUE((idty(0, 2) == scalar_t::Zero()).isFull());
0244   EXPECT_TRUE((idty(1, 2) == scalar_t::Zero()).isFull());
0245   EXPECT_TRUE((idty(2, 2) == scalar_t::One()).isFull());
0246   EXPECT_TRUE((idty(0, 3) == scalar_t::Zero()).isFull());
0247   EXPECT_TRUE((idty(1, 3) == scalar_t::Zero()).isFull());
0248   EXPECT_TRUE((idty(2, 3) == scalar_t::Zero()).isFull());
0249 
0250   // Preparatioon work
0251   vector3 z = vector::normalize(vector3{3.f, 2.f, 1.f});
0252   vector3 x = vector::normalize(vector3{2.f, -3.f, 0.f});
0253   vector3 y = vector::cross(z, x);
0254   point3 t = {2.f, 3.f, 4.f};
0255 
0256   // Test constructor from t, z, x
0257   transform3 trf1(t, z, x);
0258   ASSERT_TRUE(trf1 == trf1);
0259   transform3 trf2;
0260   trf2 = trf1;
0261 
0262   // Test printing
0263   std::cout << trf1 << std::endl;
0264 
0265   // Test comparison
0266   constexpr auto epsilon{std::numeric_limits<value_t>::epsilon()};
0267 
0268   EXPECT_TRUE(detray::algebra::approx_equal(trf1, trf1));
0269   EXPECT_TRUE(detray::algebra::approx_equal(trf1, trf1, epsilon));
0270 
0271   value_t rel_err{1.f + 10.f * epsilon};
0272   transform3 trf1_err(rel_err * t, rel_err * z, rel_err * x);
0273   EXPECT_FALSE(trf1 == trf1_err);
0274   EXPECT_TRUE(detray::algebra::approx_equal(trf1, trf1_err, 200.f * epsilon));
0275   EXPECT_FALSE(detray::algebra::approx_equal(trf1, trf1_err, 10.f * epsilon));
0276   // Cast simd vectors to different precision
0277   auto trf1_cast_f = detray::algebra::cast_to<float>(trf1);
0278   auto trf1_cast_d = detray::algebra::cast_to<double>(trf1);
0279   auto trf1_cast_i = detray::algebra::cast_to<int>(trf1);
0280 
0281   using algebra_f_t = detray::vc_soa<float>;
0282   using algebra_d_t = detray::vc_soa<double>;
0283   using algebra_i_t = detray::vc_soa<int>;
0284   static_assert(std::same_as<decltype(trf1_cast_f), dtransform3D<algebra_f_t>>);
0285   static_assert(std::same_as<decltype(trf1_cast_d), dtransform3D<algebra_d_t>>);
0286   static_assert(std::same_as<decltype(trf1_cast_i), dtransform3D<algebra_i_t>>);
0287 
0288   const auto& mat_f = trf1_cast_f.matrix();
0289   const auto& mat_d = trf1_cast_d.matrix();
0290   const auto& mat_i = trf1_cast_i.matrix();
0291   for (int j = 0; j < 3; ++j) {
0292     for (int i = 0; i < 3; ++i) {
0293       const auto& elem_ij = trf1.matrix()[i][j];
0294       EXPECT_TRUE(
0295           (mat_f[i][j] == detray::algebra::cast_to<float>(elem_ij)).isFull());
0296       EXPECT_TRUE(
0297           (mat_d[i][j] == detray::algebra::cast_to<double>(elem_ij)).isFull());
0298       EXPECT_TRUE(
0299           (mat_i[i][j] == detray::algebra::cast_to<int>(elem_ij)).isFull());
0300     }
0301   }
0302 
0303   EXPECT_TRUE((trf2(0, 0) == x[0]).isFull());
0304   EXPECT_TRUE((trf2(1, 0) == x[1]).isFull());
0305   EXPECT_TRUE((trf2(2, 0) == x[2]).isFull());
0306   EXPECT_TRUE((trf2(0, 1) == y[0]).isFull());
0307   EXPECT_TRUE((trf2(1, 1) == y[1]).isFull());
0308   EXPECT_TRUE((trf2(2, 1) == y[2]).isFull());
0309   EXPECT_TRUE((trf2(0, 2) == z[0]).isFull());
0310   EXPECT_TRUE((trf2(1, 2) == z[1]).isFull());
0311   EXPECT_TRUE((trf2(2, 2) == z[2]).isFull());
0312   EXPECT_TRUE((trf2(0, 3) == 2.f * scalar_t::One()).isFull());
0313   EXPECT_TRUE((trf2(1, 3) == 3.f * scalar_t::One()).isFull());
0314   EXPECT_TRUE((trf2(2, 3) == 4.f * scalar_t::One()).isFull());
0315 
0316   // Check that local origin translates into global translation
0317   point3 lzero = {0.f, 0.f, 0.f};
0318   point3 gzero = trf2.point_to_global(lzero);
0319   EXPECT_TRUE((gzero[0] == t[0]).isFull());
0320   EXPECT_TRUE((gzero[1] == t[1]).isFull());
0321   EXPECT_TRUE((gzero[2] == t[2]).isFull());
0322 
0323   // Check a round trip for point
0324   point3 loc_pt = {3.f, 4.f, 5.f};
0325   point3 glob_pt = trf2.point_to_global(loc_pt);
0326   point3 loc_pt_r = trf2.point_to_local(glob_pt);
0327   EXPECT_NEAR(loc_pt[0][0], loc_pt_r[0][0], tol);
0328   EXPECT_NEAR(loc_pt[1][0], loc_pt_r[1][0], tol);
0329   EXPECT_NEAR(loc_pt[2][0], loc_pt_r[2][0], tol);
0330 
0331   // Check a point versus vector transform
0332   // vector should not change if transformed by a pure translation
0333   transform3 ttrf(t);
0334 
0335   vector3 glob_vec = {1.f, 1.f, 1.f};
0336   vector3 loc_vec = ttrf.vector_to_local(glob_vec);
0337   EXPECT_NEAR(glob_vec[0][0], loc_vec[0][0], tol);
0338   EXPECT_NEAR(glob_vec[1][0], loc_vec[1][0], tol);
0339   EXPECT_NEAR(glob_vec[2][0], loc_vec[2][0], tol);
0340 
0341   // Check a round trip for vector
0342   vector3 loc_vecB = {7.f, 8.f, 9.f};
0343   vector3 glob_vecB = trf2.vector_to_local(loc_vecB);
0344   vector3 loc_vecC = trf2.vector_to_global(glob_vecB);
0345   EXPECT_NEAR(loc_vecB[0][0], loc_vecC[0][0], tol);
0346   EXPECT_NEAR(loc_vecB[1][0], loc_vecC[1][0], tol);
0347   EXPECT_NEAR(loc_vecB[2][0], loc_vecC[2][0], tol);
0348 }
0349 
0350 /// This test an SoA (Vc::Vector) based 2x3 matrix
0351 TEST(detray_algebra_soa, matrix3) {
0352   using test_algebra_t = detray::vc_soa<value_t>;
0353   using matrix_2x3_t = dmatrix<test_algebra_t, 2, 3>;
0354 
0355   // Test type traits
0356   static_assert(
0357       std::is_same_v<detray::traits::index_t<matrix_2x3_t>, std::size_t>);
0358   static_assert(std::is_same_v<detray::traits::value_t<matrix_2x3_t>, value_t>);
0359   static_assert(std::is_same_v<detray::traits::scalar_t<matrix_2x3_t>,
0360                                Vc::Vector<value_t>>);
0361   static_assert(std::is_same_v<detray::traits::vector_t<matrix_2x3_t>,
0362                                dvector2D<test_algebra_t>>);
0363 
0364   static_assert(detray::traits::rows<matrix_2x3_t> == 2);
0365   static_assert(detray::traits::columns<matrix_2x3_t> == 3);
0366   static_assert(detray::traits::max_rank<matrix_2x3_t> == 2);
0367   static_assert(detray::traits::size<matrix_2x3_t> == 6);
0368   static_assert(!detray::traits::is_square<matrix_2x3_t>);
0369   static_assert(detray::traits::is_square<dmatrix<test_algebra_t, 2, 2>>);
0370   static_assert(detray::traits::is_square<dmatrix<test_algebra_t, 3, 3>>);
0371 }
0372 
0373 /// This test an SoA (Vc::Vector) based 6x4 matrix
0374 TEST(detray_algebra_soa, matrix64) {
0375   using test_algebra_t = detray::vc_soa<value_t>;
0376 
0377   // Create the matrix.
0378   using matrix_6x4_t = dmatrix<test_algebra_t, 6, 4>;
0379   using scalar_t = dscalar<test_algebra_t>;
0380 
0381   matrix_6x4_t m;
0382 
0383   // Test type traits
0384   static_assert(
0385       std::is_same_v<detray::traits::index_t<matrix_6x4_t>, std::size_t>);
0386   static_assert(std::is_same_v<detray::traits::value_t<matrix_6x4_t>, value_t>);
0387   static_assert(std::is_same_v<detray::traits::scalar_t<matrix_6x4_t>,
0388                                Vc::Vector<value_t>>);
0389 
0390   static_assert(detray::traits::rows<matrix_6x4_t> == 6);
0391   static_assert(detray::traits::columns<matrix_6x4_t> == 4);
0392   static_assert(detray::traits::max_rank<matrix_6x4_t> == 4);
0393   static_assert(detray::traits::size<matrix_6x4_t> == 24);
0394   static_assert(!detray::traits::is_square<matrix_6x4_t>);
0395   static_assert(detray::traits::is_square<dmatrix<test_algebra_t, 4, 4>>);
0396   static_assert(detray::traits::is_square<dmatrix<test_algebra_t, 6, 6>>);
0397 
0398   // Test printing
0399   std::cout << m << std::endl;
0400 
0401   auto I64 = detray::matrix::identity<matrix_6x4_t>();
0402 
0403   // Test comparison
0404   constexpr auto epsilon{std::numeric_limits<value_t>::epsilon()};
0405 
0406   EXPECT_TRUE(detray::algebra::approx_equal(m, m));
0407   EXPECT_TRUE(detray::algebra::approx_equal(m, m, epsilon));
0408   EXPECT_FALSE(detray::algebra::approx_equal(m, I64));
0409 
0410   value_t rel_err{1.f + 10.f * epsilon};
0411   matrix_6x4_t I64_err = scalar_t(rel_err) * I64;
0412   EXPECT_FALSE(I64 == I64_err);
0413   EXPECT_TRUE(detray::algebra::approx_equal(I64, I64_err, 11.f * epsilon));
0414   EXPECT_FALSE(detray::algebra::approx_equal(I64, I64_err, 9.f * epsilon));
0415   // Cast simd vectors to different precision
0416   auto m_cast_f = detray::algebra::cast_to<float>(m);
0417   auto m_cast_d = detray::algebra::cast_to<double>(m);
0418   auto m_cast_i = detray::algebra::cast_to<int>(m);
0419 
0420   using algebra_f_t = detray::vc_soa<float>;
0421   using algebra_d_t = detray::vc_soa<double>;
0422   using algebra_i_t = detray::vc_soa<int>;
0423   static_assert(std::same_as<decltype(m_cast_f), dmatrix<algebra_f_t, 6, 4>>);
0424   static_assert(std::same_as<decltype(m_cast_d), dmatrix<algebra_d_t, 6, 4>>);
0425   static_assert(std::same_as<decltype(m_cast_i), dmatrix<algebra_i_t, 6, 4>>);
0426 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!