Acts/EventData/SubspaceHelpers.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2024 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 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/Types.hpp"
0014 #include "Acts/Utilities/AlgebraHelpers.hpp"
0015 #include "Acts/Utilities/Enumerate.hpp"
0016
0017 #include <bitset>
0018 #include <span>
0019
0020 #include <boost/container/static_vector.hpp>
0021
0022 namespace Acts {
0023
0024 /// @brief Check subspace indices for consistency
0025 ///
0026 /// Indices must be unique and within the full size of the subspace
0027 ///
0028 /// @tparam Container type of the container
0029 ///
0030 /// @param container the container of indices
0031 /// @param fullSize the full size of the subspace
0032 /// @param subspaceSize the size of the subspace
0033 ///
0034 /// @return true if the indices are consistent
0035 template <typename Container>
0036 inline static bool checkSubspaceIndices(const Container& container,
0037                                         std::size_t fullSize,
0038                                         std::size_t subspaceSize) {
0039   if (subspaceSize > fullSize) {
0040     return false;
0041   }
0042   if (container.size() != subspaceSize) {
0043     return false;
0044   }
0045   for (auto it = container.begin(); it != container.end();) {
0046     auto index = *it;
0047     if (index >= fullSize) {
0048       return false;
0049     }
0050     ++it;
0051     if (std::find(it, container.end(), index) != container.end()) {
0052       return false;
0053     }
0054   }
0055   return true;
0056 }
0057
0058 /// Serialize subspace indices to a single 64 bit integer
0059 ///
0060 /// @tparam FullSize the full size of the subspace
0061 ///
0062 /// @param indices the subspace indices
0063 ///
0064 /// @return the serialized subspace indices
0065 template <std::size_t FullSize>
0066 inline static SerializedSubspaceIndices serializeSubspaceIndices(
0067     const SubspaceIndices<FullSize>& indices)
0068   requires(FullSize <= 8)
0069 {
0070   SerializedSubspaceIndices result = 0;
0071   for (std::size_t i = 0; i < FullSize; ++i) {
0072     result |= static_cast<SerializedSubspaceIndices>(indices[i]) << (i * 8);
0073   }
0074   return result;
0075 }
0076
0077 /// Deserialize subspace indices from a single 64 bit integer
0078 ///
0079 /// @tparam FullSize the full size of the subspace
0080 ///
0081 /// @param serialized the serialized subspace indices
0082 ///
0083 /// @return the subspace indices
0084 template <std::size_t FullSize>
0085 inline static SubspaceIndices<FullSize> deserializeSubspaceIndices(
0086     SerializedSubspaceIndices serialized)
0087   requires(FullSize <= 8)
0088 {
0089   SubspaceIndices<FullSize> result;
0090   for (std::size_t i = 0; i < FullSize; ++i) {
0091     result[i] = static_cast<std::uint8_t>(serialized >> (i * 8));
0092   }
0093   return result;
0094 }
0095
0096 namespace detail {
0097
0098 /// Helper base class for subspace operations which provides common
0099 /// functionality for fixed and variable subspace helpers
0100 ///
0101 /// @tparam Derived the derived type
0102 /// @tparam FullSize the full size of the subspace
0103 template <typename Derived, std::size_t FullSize>
0104 class SubspaceHelperBase {
0105  public:
0106   static constexpr std::size_t kFullSize = FullSize;
0107
0108   using FullSquareMatrix = ActsSquareMatrix<kFullSize>;
0109
0110   bool empty() const { return self().empty(); }
0111   std::size_t size() const { return self().size(); }
0112
0113   auto operator[](std::size_t i) const { return self()[i]; }
0114
0115   auto begin() const { return self().begin(); }
0116   auto end() const { return self().end(); }
0117
0118   FullSquareMatrix fullProjector() const {
0119     FullSquareMatrix result = FullSquareMatrix::Zero();
0120     for (auto [i, index] : enumerate(*this)) {
0121       result(i, index) = 1;
0122     }
0123     return result;
0124   }
0125
0126   FullSquareMatrix fullExpander() const {
0127     FullSquareMatrix result = FullSquareMatrix::Zero();
0128     for (auto [i, index] : enumerate(*this)) {
0129       result(index, i) = 1;
0130     }
0131     return result;
0132   }
0133
0134   ProjectorBitset projectorBitset() const
0135     requires(kFullSize <= 8)
0136   {
0137     return matrixToBitset(fullProjector()).to_ullong();
0138   }
0139
0140  private:
0141   const Derived& self() const { return static_cast<const Derived&>(*this); }
0142 };
0143
0144 }  // namespace detail
0145
0146 /// Helper class for variable subspace operations
0147 ///
0148 /// @tparam FullSize the full size of the subspace
0149 /// @tparam index_t the index type
0150 template <std::size_t FullSize, typename index_t = std::uint8_t>
0151 class VariableSubspaceHelper
0152     : public detail::SubspaceHelperBase<
0153           VariableSubspaceHelper<FullSize, index_t>, FullSize> {
0154  public:
0155   static constexpr std::size_t kFullSize = FullSize;
0156
0157   using IndexType = index_t;
0158   using Container = boost::container::static_vector<IndexType, FullSize>;
0159
0160   template <typename OtherContainer>
0161   explicit VariableSubspaceHelper(const OtherContainer& indices) {
0162     assert(checkSubspaceIndices(indices, kFullSize, indices.size()) &&
0163            "Invalid indices");
0164     m_indices.resize(indices.size());
0165     std::transform(indices.begin(), indices.end(), m_indices.begin(),
0166                    [](auto index) { return static_cast<IndexType>(index); });
0167   }
0168
0169   bool empty() const { return m_indices.empty(); }
0170   std::size_t size() const { return m_indices.size(); }
0171
0172   IndexType operator[](std::size_t i) const { return m_indices[i]; }
0173
0174   auto begin() const { return m_indices.begin(); }
0175   auto end() const { return m_indices.end(); }
0176
0177  private:
0178   Container m_indices;
0179 };
0180
0181 /// Helper class for fixed subspace operations
0182 ///
0183 /// @tparam FullSize the full size of the subspace
0184 /// @tparam SubspaceSize the size of the subspace
0185 /// @tparam index_t the index type
0186 template <std::size_t FullSize, std::size_t SubspaceSize,
0187           typename index_t = std::uint8_t>
0188 class FixedSubspaceHelper
0189     : public detail::SubspaceHelperBase<
0190           FixedSubspaceHelper<FullSize, SubspaceSize, index_t>, FullSize> {
0191  public:
0192   static constexpr std::size_t kFullSize = FullSize;
0193   static constexpr std::size_t kSubspaceSize = SubspaceSize;
0194   static_assert(kSubspaceSize <= kFullSize, "Invalid subspace size");
0195
0196   using Projector = ActsMatrix<kSubspaceSize, kFullSize>;
0197   using Expander = ActsMatrix<kFullSize, kSubspaceSize>;
0198   using Vector = ActsVector<kSubspaceSize>;
0199   using SquareMatrix = ActsSquareMatrix<kSubspaceSize>;
0200   template <std::size_t K>
0201   using ApplyLeftResult = ActsMatrix<kSubspaceSize, kSubspaceSize>;
0202   template <std::size_t N>
0203   using ApplyRightResult = ActsMatrix<kSubspaceSize, kSubspaceSize>;
0204
0205   using IndexType = index_t;
0206   using Container = std::array<IndexType, kSubspaceSize>;
0207
0208   template <typename OtherContainer>
0209   explicit FixedSubspaceHelper(const OtherContainer& indices) {
0210     assert(checkSubspaceIndices(indices, kFullSize, kSubspaceSize) &&
0211            "Invalid indices");
0212     std::transform(indices.begin(), indices.end(), m_indices.begin(),
0213                    [](auto index) { return static_cast<IndexType>(index); });
0214   }
0215
0216   bool empty() const { return m_indices.empty(); }
0217   std::size_t size() const { return m_indices.size(); }
0218
0219   IndexType operator[](std::uint32_t i) const { return m_indices[i]; }
0220
0221   auto begin() const { return m_indices.begin(); }
0222   auto end() const { return m_indices.end(); }
0223
0224   Projector projector() const {
0225     Projector result = Projector::Zero();
0226     for (auto [i, index] : enumerate(*this)) {
0227       result(i, index) = 1;
0228     }
0229     return result;
0230   }
0231
0232   Expander expander() const {
0233     Expander result = Expander::Zero();
0234     for (auto [i, index] : enumerate(*this)) {
0235       result(index, i) = 1;
0236     }
0237     return result;
0238   }
0239
0240   template <std::size_t K, typename Derived>
0241   ApplyLeftResult<K> applyLeftOf(
0242       const Eigen::DenseBase<Derived>& matrix) const {
0243     assert(matrix.rows() == kFullSize && "Invalid matrix size");
0244     ApplyLeftResult<K> result = ApplyLeftResult<K>::Zero();
0245     for (auto [i, indexI] : enumerate(*this)) {
0246       for (std::size_t j = 0; j < K; ++j) {
0247         result(i, j) = matrix(indexI, j);
0248       }
0249     }
0250     return result;
0251   }
0252
0253   template <std::size_t N, typename Derived>
0254   ApplyRightResult<N> applyRightOf(
0255       const Eigen::DenseBase<Derived>& matrix) const {
0256     assert(matrix.cols() == kSubspaceSize && "Invalid matrix size");
0257     ApplyRightResult<N> result = ApplyRightResult<N>::Zero();
0258     for (std::size_t i = 0; i < N; ++i) {
0259       for (auto [j, indexJ] : enumerate(*this)) {
0260         result(i, j) = matrix(i, indexJ);
0261       }
0262     }
0263     return result;
0264   }
0265
0266   template <typename Derived>
0267   Vector projectVector(const Eigen::DenseBase<Derived>& fullVector) const {
0268     assert(fullVector.size() == kFullSize && "Invalid full vector size");
0269     Vector result = Vector::Zero();
0270     for (auto [i, index] : enumerate(*this)) {
0271       result(i) = fullVector(index);
0272     }
0273     return result;
0274   }
0275
0276   template <typename Derived>
0277   SquareMatrix projectMatrix(
0278       const Eigen::DenseBase<Derived>& fullMatrix) const {
0279     assert(fullMatrix.rows() == kFullSize && fullMatrix.cols() == kFullSize &&
0280            "Invalid full matrix size");
0281     SquareMatrix result = SquareMatrix::Zero();
0282     for (auto [i, indexI] : enumerate(*this)) {
0283       for (auto [j, indexJ] : enumerate(*this)) {
0284         result(i, j) = fullMatrix(indexI, indexJ);
0285       }
0286     }
0287     return result;
0288   }
0289
0290  private:
0291   Container m_indices{};
0292 };
0293
0294 template <std::size_t SubspaceSize>
0295 using FixedBoundSubspaceHelper =
0296     FixedSubspaceHelper<Acts::eBoundSize, SubspaceSize, std::uint8_t>;
0297 using VariableBoundSubspaceHelper =
0298     VariableSubspaceHelper<Acts::eBoundSize, std::uint8_t>;
0299
0300 /// Convert a projector to subspace indices
0301 ///
0302 /// @tparam kFullSize the full size of the subspace
0303 /// @tparam Derived the derived type
0304 ///
0305 /// @param projector the projector
0306 ///
0307 /// @return the subspace indices
0308 template <std::size_t kFullSize, typename Derived>
0309 SubspaceIndices<kFullSize> projectorToSubspaceIndices(
0310     const Eigen::DenseBase<Derived>& projector) {
0311   auto rows = static_cast<std::size_t>(projector.rows());
0312   auto cols = static_cast<std::size_t>(projector.cols());
0313   assert(cols == kFullSize && rows <= kFullSize && "Invalid projector size");
0314   SubspaceIndices<kFullSize> result;
0315   result.fill(kFullSize);
0316   for (std::size_t i = 0; i < rows; ++i) {
0317     for (std::size_t j = 0; j < cols; ++j) {
0318       assert((projector(i, j) == 0 || projector(i, j) == 1) &&
0319              "Invalid projector value");
0320       if (projector(i, j) == 1) {
0321         result[i] = j;
0322       }
0323     }
0324   }
0325   return result;
0326 }
0327
0328 }  // namespace Acts