Merge branch 'main' into dependabot/github_actions/SonarSource/sonarc…

…loud-github-c-cpp-3

Core/include/Acts/Utilities/detail/Subspace.hpp

@@ -93,7 +93,8 @@ class FixedSizeSubspace {
   /// @tparam index_t Input index type, must be convertible to std::uint8_t
   /// @param indices Unique, ordered indices
   template <typename index_t>
-  constexpr FixedSizeSubspace(const std::array<index_t, kSize>& indices) {
+  constexpr explicit FixedSizeSubspace(
+      const std::array<index_t, kSize>& indices) {
     for (std::size_t i = 0u; i < kSize; ++i) {
       assert((indices[i] < kFullSize) &&
              "Axis indices must be within the full space");
@@ -106,18 +107,18 @@ class FixedSizeSubspace {
       m_axes[i] = static_cast<std::uint8_t>(indices[i]);
     }
   }
-  // The subset can not be constructed w/o defining its axis indices.
-  FixedSizeSubspace() = delete;
-  FixedSizeSubspace(const FixedSizeSubspace&) = default;
-  FixedSizeSubspace(FixedSizeSubspace&&) = default;
-  FixedSizeSubspace& operator=(const FixedSizeSubspace&) = default;
-  FixedSizeSubspace& operator=(FixedSizeSubspace&&) = default;
 
   /// Size of the subspace.
   static constexpr std::size_t size() { return kSize; }
   /// Size of the full vector space.
   static constexpr std::size_t fullSize() { return kFullSize; }
 
+  /// Access axis index by position.
+  ///
+  /// @param i Position in the subspace
+  /// @return Axis index in the full space
+  constexpr std::size_t operator[](std::size_t i) const { return m_axes[i]; }
+
   /// Axis indices that comprise the subspace.
   ///
   /// The specific container and index type should be considered an
@@ -132,8 +133,8 @@ class FixedSizeSubspace {
     bool isContained = false;
     // always iterate over all elements to avoid branching and hope the compiler
     // can optimise this for us.
-    for (auto a : m_axes) {
-      isContained = (isContained || (a == index));
+    for (std::uint8_t a : m_axes) {
+      isContained = isContained || (a == index);
     }
     return isContained;
   }
@@ -146,12 +147,12 @@ class FixedSizeSubspace {
   ///
   /// @note Always returns a column vector regardless of the input
   template <typename fullspace_vector_t>
-  auto projectVector(const Eigen::MatrixBase<fullspace_vector_t>& full) const
-      -> SubspaceVectorFor<fullspace_vector_t> {
+  SubspaceVectorFor<fullspace_vector_t> projectVector(
+      const Eigen::MatrixBase<fullspace_vector_t>& full) const {
     EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(fullspace_vector_t, kFullSize);
 
     SubspaceVectorFor<fullspace_vector_t> sub;
-    for (auto i = 0u; i < kSize; ++i) {
+    for (std::size_t i = 0u; i < kSize; ++i) {
       sub[i] = full[m_axes[i]];
     }
     return sub;
@@ -165,13 +166,13 @@ class FixedSizeSubspace {
   ///
   /// @note Always returns a column vector regardless of the input
   template <typename subspace_vector_t>
-  auto expandVector(const Eigen::MatrixBase<subspace_vector_t>& sub) const
-      -> FullspaceVectorFor<subspace_vector_t> {
+  FullspaceVectorFor<subspace_vector_t> expandVector(
+      const Eigen::MatrixBase<subspace_vector_t>& sub) const {
     EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(subspace_vector_t, kSize);
 
     FullspaceVectorFor<subspace_vector_t> full;
     full.setZero();
-    for (auto i = 0u; i < kSize; ++i) {
+    for (std::size_t i = 0u; i < kSize; ++i) {
       full[m_axes[i]] = sub[i];
     }
     return full;
@@ -181,10 +182,10 @@ class FixedSizeSubspace {
   ///
   /// @tparam scalar_t Scalar type for the projection matrix
   template <typename scalar_t>
-  auto projector() const -> ProjectionMatrix<scalar_t> {
+  ProjectionMatrix<scalar_t> projector() const {
     ProjectionMatrix<scalar_t> proj;
     proj.setZero();
-    for (auto i = 0u; i < kSize; ++i) {
+    for (std::size_t i = 0u; i < kSize; ++i) {
       proj(i, m_axes[i]) = 1;
     }
     return proj;
@@ -194,14 +195,147 @@ class FixedSizeSubspace {
   ///
   /// @tparam scalar_t Scalar type of the generated expansion matrix
   template <typename scalar_t>
-  auto expander() const -> ExpansionMatrix<scalar_t> {
+  ExpansionMatrix<scalar_t> expander() const {
     ExpansionMatrix<scalar_t> expn;
     expn.setZero();
-    for (auto i = 0u; i < kSize; ++i) {
+    for (std::size_t i = 0u; i < kSize; ++i) {
+      expn(m_axes[i], i) = 1;
+    }
+    return expn;
+  }
+};
+
+/// Variable-size subspace representation.
+///
+/// @tparam kFullSize Size of the full vector space
+/// @tparam kSize Size of the subspace
+template <std::size_t kFullSize>
+class VariableSizeSubspace {
+  static_assert(kFullSize <= static_cast<std::size_t>(UINT8_MAX),
+                "Full vector space size is larger than the supported range");
+
+  template <typename scalar_t>
+  using FullProjectionMatrix = Eigen::Matrix<scalar_t, kFullSize, kFullSize>;
+  template <typename scalar_t>
+  using FullExpansionMatrix = Eigen::Matrix<scalar_t, kFullSize, kFullSize>;
+
+  std::size_t m_size{};
+
+  // the functionality could also be implemented using a std::bitset where each
+  // bit corresponds to an axis in the fullspace and set bits indicate which
+  // bits make up the subspace. this would be a more compact representation but
+  // complicates the implementation since we can not easily iterate over the
+  // indices of the subspace. storing the subspace indices directly requires a
+  // bit more memory but is easier to work with. for our typical use cases with
+  // n<=8, this still takes only 64bit of memory.
+  std::array<std::uint8_t, kFullSize> m_axes{};
+
+ public:
+  /// Construct from a container of axis indices.
+  ///
+  /// @tparam index_t Input index type, must be convertible to std::uint8_t
+  /// @param indices Unique, ordered indices
+  template <typename index_t, std::size_t kSize>
+  constexpr explicit VariableSizeSubspace(
+      const std::array<index_t, kSize>& indices) {
+    m_size = kSize;
+    for (std::size_t i = 0u; i < kSize; ++i) {
+      assert((indices[i] < kFullSize) &&
+             "Axis indices must be within the full space");
+      if (0u < i) {
+        assert((indices[i - 1u] < indices[i]) &&
+               "Axis indices must be unique and ordered");
+      }
+    }
+    for (std::size_t i = 0; i < kSize; ++i) {
+      m_axes[i] = static_cast<std::uint8_t>(indices[i]);
+    }
+  }
+
+  /// Size of the subspace.
+  constexpr std::size_t size() const { return m_size; }
+  /// Size of the full vector space.
+  static constexpr std::size_t fullSize() { return kFullSize; }
+
+  /// Access axis index by position.
+  ///
+  /// @param i Position in the subspace
+  /// @return Axis index in the full space
+  constexpr std::size_t operator[](std::size_t i) const {
+    assert(i < m_size);
+    return m_axes[i];
+  }
+
+  /// Check if the given axis index in the full space is part of the subspace.
+  constexpr bool contains(std::size_t index) const {
+    bool isContained = false;
+    // always iterate over all elements to avoid branching and hope the compiler
+    // can optimise this for us.
+    for (std::size_t i = 0; i < kFullSize; ++i) {
+      isContained = isContained || ((i < m_size) && (m_axes[i] == index));
+    }
+    return isContained;
+  }
+
+  /// Projection matrix that maps from the full space into the subspace.
+  ///
+  /// @tparam scalar_t Scalar type for the projection matrix
+  template <typename scalar_t>
+  FullProjectionMatrix<scalar_t> fullProjector() const {
+    FullProjectionMatrix<scalar_t> proj;
+    proj.setZero();
+    for (std::size_t i = 0u; i < m_size; ++i) {
+      proj(i, m_axes[i]) = 1;
+    }
+    return proj;
+  }
+
+  /// Expansion matrix that maps from the subspace into the full space.
+  ///
+  /// @tparam scalar_t Scalar type of the generated expansion matrix
+  template <typename scalar_t>
+  FullExpansionMatrix<scalar_t> fullExpander() const {
+    FullExpansionMatrix<scalar_t> expn;
+    expn.setZero();
+    for (std::size_t i = 0u; i < m_size; ++i) {
       expn(m_axes[i], i) = 1;
     }
     return expn;
   }
+
+  std::uint64_t projectorBits() const {
+    std::uint64_t result = 0;
+
+    for (std::size_t i = 0; i < m_size; ++i) {
+      for (std::size_t j = 0; j < kFullSize; ++j) {
+        // the bit order is defined in `Acts/Utilities/AlgebraHelpers.hpp`
+        // in `matrixToBitset`
+        std::size_t index = m_size * kFullSize - 1 - (i + j * m_size);
+        if (m_axes[i] == j) {
+          result |= (1ull << index);
+        }
+      }
+    }
+
+    return result;
+  }
+
+  std::uint64_t fullProjectorBits() const {
+    std::uint64_t result = 0;
+
+    for (std::size_t i = 0; i < kFullSize; ++i) {
+      for (std::size_t j = 0; j < kFullSize; ++j) {
+        // the bit order is defined in `Acts/Utilities/AlgebraHelpers.hpp`
+        // in `matrixToBitset`
+        std::size_t index = kFullSize * kFullSize - 1 - (i + j * kFullSize);
+        if (i < m_size && m_axes[i] == j) {
+          result |= (1ull << index);
+        }
+      }
+    }
+
+    return result;
+  }
 };
 
 /// @}

Tests/UnitTests/Core/Utilities/SubspaceTests.cpp

@@ -9,11 +9,14 @@
 #include <boost/test/unit_test.hpp>
 
 #include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/Utilities/AlgebraHelpers.hpp"
 #include "Acts/Utilities/detail/Subspace.hpp"
 
 #include <algorithm>
 #include <array>
 #include <cstddef>
+#include <cstdint>
 #include <numeric>
 #include <ostream>
 #include <tuple>
@@ -95,7 +98,7 @@ BOOST_AUTO_TEST_SUITE(UtilitiesSubspace)
 //   }
 // }
 
-BOOST_AUTO_TEST_CASE_TEMPLATE(FixedSizeSubspaceFloat, ScalarAndSubspace,
+BOOST_AUTO_TEST_CASE_TEMPLATE(FixedSizeSubspace, ScalarAndSubspace,
                               ScalarsAndFixedSizeSubspaces) {
   // extract the test types
   using Scalar = std::tuple_element_t<0, ScalarAndSubspace>;
@@ -137,4 +140,49 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(FixedSizeSubspaceFloat, ScalarAndSubspace,
   } while (std::next_permutation(fullIndices.begin(), fullIndices.end()));
 }
 
+BOOST_AUTO_TEST_CASE_TEMPLATE(VariableSizeSubspace, ScalarAndSubspace,
+                              ScalarsAndFixedSizeSubspaces) {
+  // extract the test types
+  using Scalar = std::tuple_element_t<0, ScalarAndSubspace>;
+  using FixedSubspace = std::tuple_element_t<1, ScalarAndSubspace>;
+  using VariableSubspace =
+      detail::VariableSizeSubspace<FixedSubspace::fullSize()>;
+
+  auto fullIndices = makeMonotonicIndices(FixedSubspace::fullSize());
+
+  // in principle, we would like to iterate over all possible ordered subsets
+  // from the full space indices with a size identical to the subspace. since i
+  // do not know how to do that in a simple manner, we are iterating over all
+  // permutations of the full indices and pick the first n elements. this should
+  // give a reasonable set of different subspace configurations.
+  do {
+    auto indices = selectFixedIndices<FixedSubspace::size()>(fullIndices);
+    FixedSubspace fixedSubspace(indices);
+    VariableSubspace variableSubspace(indices);
+
+    BOOST_CHECK_EQUAL(variableSubspace.size(), fixedSubspace.size());
+    BOOST_CHECK_EQUAL(variableSubspace.fullSize(), fixedSubspace.fullSize());
+
+    auto fixedProjector = fixedSubspace.template projector<Scalar>();
+    std::uint64_t fixedProjectorBits =
+        matrixToBitset(fixedProjector).to_ullong();
+
+    Eigen::Matrix<Scalar, FixedSubspace::fullSize(), FixedSubspace::fullSize()>
+        fixedFullProjector;
+    fixedFullProjector.setZero();
+    fixedFullProjector.template topLeftCorner<FixedSubspace::size(),
+                                              FixedSubspace::fullSize()>() =
+        fixedProjector;
+    std::uint64_t fixedFullProjectorBits =
+        matrixToBitset(fixedFullProjector).to_ullong();
+
+    std::uint64_t variableProjectorBits = variableSubspace.projectorBits();
+    std::uint64_t variableFullProjectorBits =
+        variableSubspace.fullProjectorBits();
+
+    BOOST_CHECK_EQUAL(variableProjectorBits, fixedProjectorBits);
+    BOOST_CHECK_EQUAL(variableFullProjectorBits, fixedFullProjectorBits);
+  } while (std::next_permutation(fullIndices.begin(), fullIndices.end()));
+}
+
 BOOST_AUTO_TEST_SUITE_END()