Acts/EventData/SubspaceHelpers.hpp

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 #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 <ranges>
0018
0019 #include <boost/container/static_vector.hpp>
0020
0021 namespace Acts {
0022
0023 /// @brief Check subspace indices for consistency
0024 ///
0025 /// Indices must be unique and within the full size of the subspace
0026 ///
0027 /// @tparam index_range_t the type of the container of indices
0028 ///
0029 /// @param indexRange the range of indices
0030 /// @param fullSize the full size of the subspace
0031 /// @param subspaceSize the size of the subspace
0032 ///
0033 /// @return true if the indices are consistent
0034 template <std::ranges::sized_range index_range_t>
0035 inline static bool checkSubspaceIndices(const index_range_t& indexRange,
0036                                         std::size_t fullSize,
0037                                         std::size_t subspaceSize) {
0038   if (subspaceSize > fullSize) {
0039     return false;
0040   }
0041   if (static_cast<std::size_t>(indexRange.size()) != subspaceSize) {
0042     return false;
0043   }
0044   for (auto it = indexRange.begin(); it != indexRange.end();) {
0045     auto index = *it;
0046     if (index >= fullSize) {
0047       return false;
0048     }
0049     ++it;
0050     if (std::find(it, indexRange.end(), index) != indexRange.end()) {
0051       return false;
0052     }
0053   }
0054   return true;
0055 }
0056
0057 /// Serialize subspace indices to a single 64 bit integer
0058 ///
0059 /// @tparam FullSize the full size of the subspace
0060 ///
0061 /// @param indices the subspace indices
0062 ///
0063 /// @return the serialized subspace indices
0064 template <std::size_t FullSize>
0065 inline static SerializedSubspaceIndices serializeSubspaceIndices(
0066     const SubspaceIndices<FullSize>& indices)
0067   requires(FullSize <= 8)
0068 {
0069   SerializedSubspaceIndices result = 0;
0070   for (std::size_t i = 0; i < FullSize; ++i) {
0071     result |= static_cast<SerializedSubspaceIndices>(indices[i] & 0xFF)
0072               << (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)) & 0xFF);
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   bool contains(std::uint8_t index) const {
0119     return std::find(begin(), end(), index) != end();
0120   }
0121   std::size_t indexOf(std::uint8_t index) const {
0122     auto it = std::find(begin(), end(), index);
0123     return it != end() ? std::distance(begin(), it) : kFullSize;
0124   }
0125
0126   template <typename EigenDerived>
0127   ActsVector<kFullSize> expandVector(
0128       const Eigen::DenseBase<EigenDerived>& vector) const {
0129     ActsVector<kFullSize> result = ActsVector<kFullSize>::Zero();
0130     for (auto [i, index] : enumerate(*this)) {
0131       result(index) = vector(i);
0132     }
0133     return result;
0134   }
0135
0136   template <typename EigenDerived>
0137   FullSquareMatrix expandMatrix(
0138       const Eigen::DenseBase<EigenDerived>& matrix) const {
0139     FullSquareMatrix result = FullSquareMatrix::Zero();
0140     for (auto [i, indexI] : enumerate(*this)) {
0141       for (auto [j, indexJ] : enumerate(*this)) {
0142         result(indexI, indexJ) = matrix(i, j);
0143       }
0144     }
0145     return result;
0146   }
0147
0148   FullSquareMatrix fullProjector() const {
0149     FullSquareMatrix result = FullSquareMatrix::Zero();
0150     for (auto [i, index] : enumerate(*this)) {
0151       result(i, index) = 1;
0152     }
0153     return result;
0154   }
0155
0156   FullSquareMatrix fullExpander() const {
0157     FullSquareMatrix result = FullSquareMatrix::Zero();
0158     for (auto [i, index] : enumerate(*this)) {
0159       result(index, i) = 1;
0160     }
0161     return result;
0162   }
0163
0164   ProjectorBitset projectorBitset() const
0165     requires(kFullSize <= 8)
0166   {
0167     return matrixToBitset(fullProjector()).to_ullong();
0168   }
0169
0170  private:
0171   const Derived& self() const { return static_cast<const Derived&>(*this); }
0172 };
0173
0174 }  // namespace detail
0175
0176 /// Helper class for variable subspace operations
0177 ///
0178 /// @tparam FullSize the full size of the subspace
0179 /// @tparam index_t the index type
0180 template <std::size_t FullSize, typename index_t = std::uint8_t>
0181 class VariableSubspaceHelper
0182     : public detail::SubspaceHelperBase<
0183           VariableSubspaceHelper<FullSize, index_t>, FullSize> {
0184  public:
0185   static constexpr std::size_t kFullSize = FullSize;
0186
0187   using IndexType = index_t;
0188   using Container = boost::container::static_vector<IndexType, FullSize>;
0189
0190   template <std::ranges::sized_range other_index_range_t>
0191   explicit VariableSubspaceHelper(const other_index_range_t& indices) {
0192     assert(checkSubspaceIndices(indices, kFullSize, indices.size()) &&
0193            "Invalid indices");
0194     m_indices.resize(indices.size());
0195     std::transform(indices.begin(), indices.end(), m_indices.begin(),
0196                    [](auto index) { return static_cast<IndexType>(index); });
0197   }
0198
0199   bool empty() const { return m_indices.empty(); }
0200   std::size_t size() const { return m_indices.size(); }
0201   const Container& indices() const { return m_indices; }
0202
0203   IndexType operator[](std::size_t i) const { return m_indices[i]; }
0204
0205   auto begin() const { return m_indices.begin(); }
0206   auto end() const { return m_indices.end(); }
0207
0208  private:
0209   Container m_indices;
0210 };
0211
0212 /// Helper class for fixed subspace operations
0213 ///
0214 /// @tparam FullSize the full size of the subspace
0215 /// @tparam SubspaceSize the size of the subspace
0216 /// @tparam index_t the index type
0217 template <std::size_t FullSize, std::size_t SubspaceSize,
0218           typename index_t = std::uint8_t>
0219 class FixedSubspaceHelper
0220     : public detail::SubspaceHelperBase<
0221           FixedSubspaceHelper<FullSize, SubspaceSize, index_t>, FullSize> {
0222  public:
0223   static constexpr std::size_t kFullSize = FullSize;
0224   static constexpr std::size_t kSubspaceSize = SubspaceSize;
0225   static_assert(kSubspaceSize <= kFullSize, "Invalid subspace size");
0226
0227   using Projector = ActsMatrix<kSubspaceSize, kFullSize>;
0228   using Expander = ActsMatrix<kFullSize, kSubspaceSize>;
0229   using Vector = ActsVector<kSubspaceSize>;
0230   using SquareMatrix = ActsSquareMatrix<kSubspaceSize>;
0231   template <std::size_t K>
0232   using ApplyLeftResult = ActsMatrix<kSubspaceSize, kSubspaceSize>;
0233   template <std::size_t N>
0234   using ApplyRightResult = ActsMatrix<kSubspaceSize, kSubspaceSize>;
0235
0236   using IndexType = index_t;
0237   using Container = std::array<IndexType, kSubspaceSize>;
0238
0239   template <std::ranges::sized_range other_index_range_t>
0240   explicit FixedSubspaceHelper(const other_index_range_t& indices) {
0241     assert(checkSubspaceIndices(indices, kFullSize, kSubspaceSize) &&
0242            "Invalid indices");
0243     std::transform(indices.begin(), indices.end(), m_indices.begin(),
0244                    [](auto index) { return static_cast<IndexType>(index); });
0245   }
0246
0247   bool empty() const { return m_indices.empty(); }
0248   std::size_t size() const { return m_indices.size(); }
0249   const Container& indices() const { return m_indices; }
0250
0251   IndexType operator[](std::uint32_t i) const { return m_indices[i]; }
0252
0253   auto begin() const { return m_indices.begin(); }
0254   auto end() const { return m_indices.end(); }
0255
0256   Projector projector() const {
0257     Projector result = Projector::Zero();
0258     for (auto [i, index] : enumerate(*this)) {
0259       result(i, index) = 1;
0260     }
0261     return result;
0262   }
0263
0264   Expander expander() const {
0265     Expander result = Expander::Zero();
0266     for (auto [i, index] : enumerate(*this)) {
0267       result(index, i) = 1;
0268     }
0269     return result;
0270   }
0271
0272   template <std::size_t K, typename Derived>
0273   ApplyLeftResult<K> applyLeftOf(
0274       const Eigen::DenseBase<Derived>& matrix) const {
0275     assert(matrix.rows() == kFullSize && "Invalid matrix size");
0276     ApplyLeftResult<K> result = ApplyLeftResult<K>::Zero();
0277     for (auto [i, indexI] : enumerate(*this)) {
0278       for (std::size_t j = 0; j < K; ++j) {
0279         result(i, j) = matrix(indexI, j);
0280       }
0281     }
0282     return result;
0283   }
0284
0285   template <std::size_t N, typename Derived>
0286   ApplyRightResult<N> applyRightOf(
0287       const Eigen::DenseBase<Derived>& matrix) const {
0288     assert(matrix.cols() == kSubspaceSize && "Invalid matrix size");
0289     ApplyRightResult<N> result = ApplyRightResult<N>::Zero();
0290     for (std::size_t i = 0; i < N; ++i) {
0291       for (auto [j, indexJ] : enumerate(*this)) {
0292         result(i, j) = matrix(i, indexJ);
0293       }
0294     }
0295     return result;
0296   }
0297
0298   template <typename Derived>
0299   Vector projectVector(const Eigen::DenseBase<Derived>& fullVector) const {
0300     assert(fullVector.size() == kFullSize && "Invalid full vector size");
0301     Vector result = Vector::Zero();
0302     for (auto [i, index] : enumerate(*this)) {
0303       result(i) = fullVector(index);
0304     }
0305     return result;
0306   }
0307
0308   template <typename Derived>
0309   SquareMatrix projectMatrix(
0310       const Eigen::DenseBase<Derived>& fullMatrix) const {
0311     assert(fullMatrix.rows() == kFullSize && fullMatrix.cols() == kFullSize &&
0312            "Invalid full matrix size");
0313     SquareMatrix result = SquareMatrix::Zero();
0314     for (auto [i, indexI] : enumerate(*this)) {
0315       for (auto [j, indexJ] : enumerate(*this)) {
0316         result(i, j) = fullMatrix(indexI, indexJ);
0317       }
0318     }
0319     return result;
0320   }
0321
0322  private:
0323   Container m_indices{};
0324 };
0325
0326 template <std::size_t SubspaceSize>
0327 using FixedBoundSubspaceHelper =
0328     FixedSubspaceHelper<Acts::eBoundSize, SubspaceSize, std::uint8_t>;
0329 using VariableBoundSubspaceHelper =
0330     VariableSubspaceHelper<Acts::eBoundSize, std::uint8_t>;
0331
0332 /// Convert a projector to subspace indices
0333 ///
0334 /// @tparam kFullSize the full size of the subspace
0335 /// @tparam Derived the derived type
0336 ///
0337 /// @param projector the projector
0338 ///
0339 /// @return the subspace indices
0340 template <std::size_t kFullSize, typename Derived>
0341 SubspaceIndices<kFullSize> projectorToSubspaceIndices(
0342     const Eigen::DenseBase<Derived>& projector) {
0343   auto rows = static_cast<std::size_t>(projector.rows());
0344   auto cols = static_cast<std::size_t>(projector.cols());
0345   assert(cols == kFullSize && rows <= kFullSize && "Invalid projector size");
0346   SubspaceIndices<kFullSize> result;
0347   result.fill(kFullSize);
0348   for (std::size_t i = 0; i < rows; ++i) {
0349     for (std::size_t j = 0; j < cols; ++j) {
0350       assert((projector(i, j) == 0 || projector(i, j) == 1) &&
0351              "Invalid projector value");
0352       if (projector(i, j) == 1) {
0353         result[i] = j;
0354       }
0355     }
0356   }
0357   return result;
0358 }
0359
0360 }  // namespace Acts