Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-06 08:07:28

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020 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 http://mozilla.org/MPL/2.0/.
0008 
0009 #pragma once
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 
0013 #include <array>
0014 #include <cassert>
0015 #include <cstddef>
0016 #include <cstdint>
0017 
0018 namespace Acts::detail {
0019 
0020 /// @defgroup subspace Linear subspace definitions
0021 ///
0022 /// All types in this group define a linear subspace of a larger vector space
0023 /// with the following properties: the subspace is defined by a set of axis
0024 /// indices in the full space, i.e. there is no need for a non-trivial
0025 /// projection matrix, all subspace axis indices are unique and there are no
0026 /// duplicated ones, and the set of axis indices must be ordered, i.e. the axis
0027 /// order in the subspace is the same as in the full space. The last requirement
0028 /// is not strictly required by the implementation, but is still added to
0029 /// simplify reasoning.
0030 ///
0031 /// Only the size of the subspace are defined by the types at compile time.
0032 /// Which axes comprise the subspace is defined at runtime. This allows to use
0033 /// fixed-size computations (selected at compile time) but avoids the
0034 /// combinatorial explosion of also having to handle all possible combination of
0035 /// axes at compile time. This was tried previously and resulted in sizable
0036 /// resource consumption at compile time without runtime benefits.
0037 ///
0038 /// All types are intentionally using `std::size_t` as their template values,
0039 /// instead of the more specific index enums, to reduce the number of templates.
0040 /// This is fully compatible as the index enums are required to be convertible
0041 /// to `std::size_t`.
0042 ///
0043 /// All types intentionally only define the subspace but not how vectors
0044 /// and matrices are stored to avoid unnecessary coupling between components,
0045 /// i.e here between the pure definition and the storage.
0046 ///
0047 /// All types provide `.projectVector(...)` and `.exandVector(...)` methods to
0048 /// convert to/from the subspace. They also provide `.projector()` and
0049 /// `.expander()` methods to create projection and expansion matrices if they
0050 /// are required explicitly. For the specific subspace requirements listed
0051 /// above, the projection and expansion matrices are transpose to each other. In
0052 /// the general case, this does not have to be the case and users are encouraged
0053 /// to use `.projector()` and `.expander()` instead of e.g.
0054 /// `.projector().transpose()`.
0055 ///
0056 /// @{
0057 
0058 /// Fixed-size subspace representation.
0059 ///
0060 /// @tparam kFullSize Size of the full vector space
0061 /// @tparam kSize Size of the subspace
0062 template <std::size_t kFullSize, std::size_t kSize>
0063 class FixedSizeSubspace {
0064   static_assert(kFullSize <= static_cast<std::size_t>(
0065                                  std::numeric_limits<std::uint8_t>::max()),
0066                 "Full vector space size is larger than the supported range");
0067   static_assert(1u <= kSize, "Subspace size must be at least 1");
0068   static_assert(kSize <= kFullSize,
0069                 "Subspace can only be as large as the full space");
0070 
0071   template <typename source_t>
0072   using SubspaceVectorFor = Eigen::Matrix<typename source_t::Scalar, kSize, 1>;
0073   template <typename source_t>
0074   using FullspaceVectorFor =
0075       Eigen::Matrix<typename source_t::Scalar, kFullSize, 1>;
0076   template <typename scalar_t>
0077   using ProjectionMatrix = Eigen::Matrix<scalar_t, kSize, kFullSize>;
0078   template <typename scalar_t>
0079   using ExpansionMatrix = Eigen::Matrix<scalar_t, kFullSize, kSize>;
0080 
0081   // the functionality could also be implemented using a std::bitset where each
0082   // bit corresponds to an axis in the fullspace and set bits indicate which
0083   // bits make up the subspace. this would be a more compact representation but
0084   // complicates the implementation since we can not easily iterate over the
0085   // indices of the subspace. storing the subspace indices directly requires a
0086   // bit more memory but is easier to work with. for our typical use cases with
0087   // n<=8, this still takes only 64bit of memory.
0088   std::array<std::uint8_t, kSize> m_axes;
0089 
0090  public:
0091   /// Construct from a container of axis indices.
0092   ///
0093   /// @tparam index_t Input index type, must be convertible to std::uint8_t
0094   /// @param indices Unique, ordered indices
0095   template <typename index_t>
0096   constexpr explicit FixedSizeSubspace(
0097       const std::array<index_t, kSize>& indices) {
0098     for (std::size_t i = 0u; i < kSize; ++i) {
0099       assert((indices[i] < kFullSize) &&
0100              "Axis indices must be within the full space");
0101       if (0u < i) {
0102         assert((indices[i - 1u] < indices[i]) &&
0103                "Axis indices must be unique and ordered");
0104       }
0105     }
0106     for (std::size_t i = 0; i < kSize; ++i) {
0107       m_axes[i] = static_cast<std::uint8_t>(indices[i]);
0108     }
0109   }
0110 
0111   /// Size of the subspace.
0112   static constexpr std::size_t size() { return kSize; }
0113   /// Size of the full vector space.
0114   static constexpr std::size_t fullSize() { return kFullSize; }
0115 
0116   /// Access axis index by position.
0117   ///
0118   /// @param i Position in the subspace
0119   /// @return Axis index in the full space
0120   constexpr std::size_t operator[](std::size_t i) const { return m_axes[i]; }
0121 
0122   std::size_t indexOf(std::size_t axis) const {
0123     for (std::size_t i = 0; i < kSize; ++i) {
0124       if (m_axes[i] == axis) {
0125         return i;
0126       }
0127     }
0128     return kSize;
0129   }
0130 
0131   /// Axis indices that comprise the subspace.
0132   ///
0133   /// The specific container and index type should be considered an
0134   /// implementation detail. Users should treat the return type as a generic
0135   /// container whose elements are convertible to `std::size_t`.
0136   constexpr const std::array<std::uint8_t, kSize>& indices() const {
0137     return m_axes;
0138   }
0139 
0140   /// Check if the given axis index in the full space is part of the subspace.
0141   constexpr bool contains(std::size_t index) const {
0142     bool isContained = false;
0143     // always iterate over all elements to avoid branching and hope the compiler
0144     // can optimise this for us.
0145     for (std::uint8_t a : m_axes) {
0146       isContained = isContained || (a == index);
0147     }
0148     return isContained;
0149   }
0150 
0151   /// Project a full vector into the subspace.
0152   ///
0153   /// @tparam fullspace_vector_t Vector type in the full space
0154   /// @param full Vector in the full space
0155   /// @return Subspace vector w/ just the configured axis components
0156   ///
0157   /// @note Always returns a column vector regardless of the input
0158   template <typename fullspace_vector_t>
0159   SubspaceVectorFor<fullspace_vector_t> projectVector(
0160       const Eigen::MatrixBase<fullspace_vector_t>& full) const {
0161     EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(fullspace_vector_t, kFullSize);
0162 
0163     SubspaceVectorFor<fullspace_vector_t> sub;
0164     for (std::size_t i = 0u; i < kSize; ++i) {
0165       sub[i] = full[m_axes[i]];
0166     }
0167     return sub;
0168   }
0169 
0170   /// Expand a subspace vector into the full space.
0171   ///
0172   /// @tparam subspace_vector_t Subspace vector type
0173   /// @param sub Subspace vector
0174   /// @return Vector in the full space w/ zeros for non-subspace components
0175   ///
0176   /// @note Always returns a column vector regardless of the input
0177   template <typename subspace_vector_t>
0178   FullspaceVectorFor<subspace_vector_t> expandVector(
0179       const Eigen::MatrixBase<subspace_vector_t>& sub) const {
0180     EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(subspace_vector_t, kSize);
0181 
0182     FullspaceVectorFor<subspace_vector_t> full;
0183     full.setZero();
0184     for (std::size_t i = 0u; i < kSize; ++i) {
0185       full[m_axes[i]] = sub[i];
0186     }
0187     return full;
0188   }
0189 
0190   /// Projection matrix that maps from the full space into the subspace.
0191   ///
0192   /// @tparam scalar_t Scalar type for the projection matrix
0193   template <typename scalar_t>
0194   ProjectionMatrix<scalar_t> projector() const {
0195     ProjectionMatrix<scalar_t> proj;
0196     proj.setZero();
0197     for (std::size_t i = 0u; i < kSize; ++i) {
0198       proj(i, m_axes[i]) = 1;
0199     }
0200     return proj;
0201   }
0202 
0203   /// Expansion matrix that maps from the subspace into the full space.
0204   ///
0205   /// @tparam scalar_t Scalar type of the generated expansion matrix
0206   template <typename scalar_t>
0207   ExpansionMatrix<scalar_t> expander() const {
0208     ExpansionMatrix<scalar_t> expn;
0209     expn.setZero();
0210     for (std::size_t i = 0u; i < kSize; ++i) {
0211       expn(m_axes[i], i) = 1;
0212     }
0213     return expn;
0214   }
0215 };
0216 
0217 /// Variable-size subspace representation.
0218 ///
0219 /// @tparam kFullSize Size of the full vector space
0220 /// @tparam kSize Size of the subspace
0221 template <std::size_t kFullSize>
0222 class VariableSizeSubspace {
0223   static_assert(kFullSize <= static_cast<std::size_t>(UINT8_MAX),
0224                 "Full vector space size is larger than the supported range");
0225 
0226   template <typename scalar_t>
0227   using FullProjectionMatrix = Eigen::Matrix<scalar_t, kFullSize, kFullSize>;
0228   template <typename scalar_t>
0229   using FullExpansionMatrix = Eigen::Matrix<scalar_t, kFullSize, kFullSize>;
0230 
0231   std::size_t m_size{};
0232 
0233   // the functionality could also be implemented using a std::bitset where each
0234   // bit corresponds to an axis in the fullspace and set bits indicate which
0235   // bits make up the subspace. this would be a more compact representation but
0236   // complicates the implementation since we can not easily iterate over the
0237   // indices of the subspace. storing the subspace indices directly requires a
0238   // bit more memory but is easier to work with. for our typical use cases with
0239   // n<=8, this still takes only 64bit of memory.
0240   std::array<std::uint8_t, kFullSize> m_axes{};
0241 
0242  public:
0243   /// Construct from a container of axis indices.
0244   ///
0245   /// @tparam index_t Input index type, must be convertible to std::uint8_t
0246   /// @param indices Unique, ordered indices
0247   template <typename index_t, std::size_t kSize>
0248   constexpr explicit VariableSizeSubspace(
0249       const std::array<index_t, kSize>& indices) {
0250     m_size = kSize;
0251     for (std::size_t i = 0u; i < kSize; ++i) {
0252       assert((indices[i] < kFullSize) &&
0253              "Axis indices must be within the full space");
0254       if (0u < i) {
0255         assert((indices[i - 1u] < indices[i]) &&
0256                "Axis indices must be unique and ordered");
0257       }
0258     }
0259     for (std::size_t i = 0; i < kSize; ++i) {
0260       m_axes[i] = static_cast<std::uint8_t>(indices[i]);
0261     }
0262   }
0263 
0264   /// Size of the subspace.
0265   constexpr std::size_t size() const { return m_size; }
0266   /// Size of the full vector space.
0267   static constexpr std::size_t fullSize() { return kFullSize; }
0268 
0269   /// Access axis index by position.
0270   ///
0271   /// @param i Position in the subspace
0272   /// @return Axis index in the full space
0273   constexpr std::size_t operator[](std::size_t i) const {
0274     assert(i < m_size);
0275     return m_axes[i];
0276   }
0277 
0278   std::size_t indexOf(std::size_t axis) const {
0279     for (std::size_t i = 0; i < m_size; ++i) {
0280       if (m_axes[i] == axis) {
0281         return i;
0282       }
0283     }
0284     return m_size;
0285   }
0286 
0287   /// Check if the given axis index in the full space is part of the subspace.
0288   constexpr bool contains(std::size_t index) const {
0289     bool isContained = false;
0290     // always iterate over all elements to avoid branching and hope the compiler
0291     // can optimise this for us.
0292     for (std::size_t i = 0; i < kFullSize; ++i) {
0293       isContained = isContained || ((i < m_size) && (m_axes[i] == index));
0294     }
0295     return isContained;
0296   }
0297 
0298   /// Projection matrix that maps from the full space into the subspace.
0299   ///
0300   /// @tparam scalar_t Scalar type for the projection matrix
0301   template <typename scalar_t>
0302   FullProjectionMatrix<scalar_t> fullProjector() const {
0303     FullProjectionMatrix<scalar_t> proj;
0304     proj.setZero();
0305     for (std::size_t i = 0u; i < m_size; ++i) {
0306       proj(i, m_axes[i]) = 1;
0307     }
0308     return proj;
0309   }
0310 
0311   /// Expansion matrix that maps from the subspace into the full space.
0312   ///
0313   /// @tparam scalar_t Scalar type of the generated expansion matrix
0314   template <typename scalar_t>
0315   FullExpansionMatrix<scalar_t> fullExpander() const {
0316     FullExpansionMatrix<scalar_t> expn;
0317     expn.setZero();
0318     for (std::size_t i = 0u; i < m_size; ++i) {
0319       expn(m_axes[i], i) = 1;
0320     }
0321     return expn;
0322   }
0323 };
0324 
0325 /// @}
0326 
0327 }  // namespace Acts::detail