Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:27:55

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>(UINT8_MAX),
0065                 "Full vector space size is larger than the supported range");
0066   static_assert(1u <= kSize, "Subspace size must be at least 1");
0067   static_assert(kSize <= kFullSize,
0068                 "Subspace can only be as large as the full space");
0069 
0070   template <typename source_t>
0071   using SubspaceVectorFor = Eigen::Matrix<typename source_t::Scalar, kSize, 1>;
0072   template <typename source_t>
0073   using FullspaceVectorFor =
0074       Eigen::Matrix<typename source_t::Scalar, kFullSize, 1>;
0075   template <typename scalar_t>
0076   using ProjectionMatrix = Eigen::Matrix<scalar_t, kSize, kFullSize>;
0077   template <typename scalar_t>
0078   using ExpansionMatrix = Eigen::Matrix<scalar_t, kFullSize, kSize>;
0079 
0080   // the functionality could also be implemented using a std::bitset where each
0081   // bit corresponds to an axis in the fullspace and set bits indicate which
0082   // bits make up the subspace. this would be a more compact representation but
0083   // complicates the implementation since we can not easily iterate over the
0084   // indices of the subspace. storing the subspace indices directly requires a
0085   // bit more memory but is easier to work with. for our typical use cases with
0086   // n<=8, this still takes only 64bit of memory.
0087   std::array<uint8_t, kSize> m_axes;
0088 
0089  public:
0090   /// Construct from a container of axis indices.
0091   ///
0092   /// @tparam index_t Input index type, must be convertible to uint8_t
0093   /// @param indices Unique, ordered indices
0094   template <typename index_t>
0095   constexpr FixedSizeSubspace(const std::array<index_t, kSize>& indices) {
0096     for (std::size_t i = 0u; i < kSize; ++i) {
0097       assert((indices[i] < kFullSize) &&
0098              "Axis indices must be within the full space");
0099       if (0u < i) {
0100         assert((indices[i - 1u] < indices[i]) &&
0101                "Axis indices must be unique and ordered");
0102       }
0103     }
0104     for (std::size_t i = 0; i < kSize; ++i) {
0105       m_axes[i] = static_cast<uint8_t>(indices[i]);
0106     }
0107   }
0108   // The subset can not be constructed w/o defining its axis indices.
0109   FixedSizeSubspace() = delete;
0110   FixedSizeSubspace(const FixedSizeSubspace&) = default;
0111   FixedSizeSubspace(FixedSizeSubspace&&) = default;
0112   FixedSizeSubspace& operator=(const FixedSizeSubspace&) = default;
0113   FixedSizeSubspace& operator=(FixedSizeSubspace&&) = default;
0114 
0115   /// Size of the subspace.
0116   static constexpr std::size_t size() { return kSize; }
0117   /// Size of the full vector space.
0118   static constexpr std::size_t fullSize() { return kFullSize; }
0119 
0120   /// Axis indices that comprise the subspace.
0121   ///
0122   /// The specific container and index type should be considered an
0123   /// implementation detail. Users should treat the return type as a generic
0124   /// container whose elements are convertible to `std::size_t`.
0125   constexpr const std::array<uint8_t, kSize>& indices() const { return m_axes; }
0126 
0127   /// Check if the given axis index in the full space is part of the subspace.
0128   constexpr bool contains(std::size_t index) const {
0129     bool isContained = false;
0130     // always iterate over all elements to avoid branching and hope the compiler
0131     // can optimise this for us.
0132     for (auto a : m_axes) {
0133       isContained = (isContained || (a == index));
0134     }
0135     return isContained;
0136   }
0137 
0138   /// Project a full vector into the subspace.
0139   ///
0140   /// @tparam fullspace_vector_t Vector type in the full space
0141   /// @param full Vector in the full space
0142   /// @return Subspace vector w/ just the configured axis components
0143   ///
0144   /// @note Always returns a column vector regardless of the input
0145   template <typename fullspace_vector_t>
0146   auto projectVector(const Eigen::MatrixBase<fullspace_vector_t>& full) const
0147       -> SubspaceVectorFor<fullspace_vector_t> {
0148     EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(fullspace_vector_t, kFullSize);
0149 
0150     SubspaceVectorFor<fullspace_vector_t> sub;
0151     for (auto i = 0u; i < kSize; ++i) {
0152       sub[i] = full[m_axes[i]];
0153     }
0154     return sub;
0155   }
0156 
0157   /// Expand a subspace vector into the full space.
0158   ///
0159   /// @tparam subspace_vector_t Subspace vector type
0160   /// @param sub Subspace vector
0161   /// @return Vector in the full space w/ zeros for non-subspace components
0162   ///
0163   /// @note Always returns a column vector regardless of the input
0164   template <typename subspace_vector_t>
0165   auto expandVector(const Eigen::MatrixBase<subspace_vector_t>& sub) const
0166       -> FullspaceVectorFor<subspace_vector_t> {
0167     EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(subspace_vector_t, kSize);
0168 
0169     FullspaceVectorFor<subspace_vector_t> full;
0170     full.setZero();
0171     for (auto i = 0u; i < kSize; ++i) {
0172       full[m_axes[i]] = sub[i];
0173     }
0174     return full;
0175   }
0176 
0177   /// Projection matrix that maps from the full space into the subspace.
0178   ///
0179   /// @tparam scalar_t Scalar type for the projection matrix
0180   template <typename scalar_t>
0181   auto projector() const -> ProjectionMatrix<scalar_t> {
0182     ProjectionMatrix<scalar_t> proj;
0183     proj.setZero();
0184     for (auto i = 0u; i < kSize; ++i) {
0185       proj(i, m_axes[i]) = 1;
0186     }
0187     return proj;
0188   }
0189 
0190   /// Expansion matrix that maps from the subspace into the full space.
0191   ///
0192   /// @tparam scalar_t Scalar type of the generated expansion matrix
0193   template <typename scalar_t>
0194   auto expander() const -> ExpansionMatrix<scalar_t> {
0195     ExpansionMatrix<scalar_t> expn;
0196     expn.setZero();
0197     for (auto i = 0u; i < kSize; ++i) {
0198       expn(m_axes[i], i) = 1;
0199     }
0200     return expn;
0201   }
0202 };
0203 
0204 /// @}
0205 
0206 }  // namespace Acts::detail