refactor: Rework projector

Core/include/Acts/Definitions/Alignment.hpp

@@ -19,7 +19,7 @@ namespace Acts {
 /// This must be a regular `enum` and not a scoped `enum class` to allow
 /// implicit conversion to an integer. The enum value are thus visible directly
 /// in `namespace Acts` and are prefixed to avoid naming collisions.
-enum AlignmentIndices : unsigned int {
+enum AlignmentIndices : std::uint8_t {
   // Center of geometry object in global 3D cartesian coordinates
   eAlignmentCenter0 = 0u,
   eAlignmentCenter1 = eAlignmentCenter0 + 1u,

Core/include/Acts/Definitions/Common.hpp

@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include <cstdint>
 #include <iosfwd>
 
 namespace Acts {
@@ -40,7 +41,7 @@ enum NoiseUpdateMode : int { removeNoise = -1, addNoise = 1 };
 /// This index enum is not user-configurable (in contrast e.g. to the track
 /// parameter index enums) since it must be compatible with varying
 /// dimensionality (2d-4d) and other access methods (`.{x,y,z}()` accessors).
-enum CoordinateIndices : unsigned int {
+enum CoordinateIndices : std::uint8_t {
   // generic position-like access
   ePos0 = 0,
   ePos1 = 1,

Core/include/Acts/Definitions/TrackParametrization.hpp

@@ -10,6 +10,7 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 
+#include <cstdint>
 #include <type_traits>
 
 // The user can override the track parameters ordering. If the preprocessor
@@ -34,7 +35,7 @@ namespace Acts {
 /// This must be a regular `enum` and not a scoped `enum class` to allow
 /// implicit conversion to an integer. The enum value are thus visible directly
 /// in `namespace Acts` and are prefixed to avoid naming collisions.
-enum BoundIndices : unsigned int {
+enum BoundIndices : std::uint8_t {
   // Local position on the reference surface.
   // This is intentionally named different from the position components in
   // the other data vectors, to clarify that this is defined on a surface
@@ -61,7 +62,7 @@ enum BoundIndices : unsigned int {
 /// This must be a regular `enum` and not a scoped `enum class` to allow
 /// implicit conversion to an integer. The enum value are thus visible directly
 /// in `namespace Acts` and are prefixed to avoid naming collisions.
-enum FreeIndices : unsigned int {
+enum FreeIndices : std::uint8_t {
   // Spatial position
   // The spatial position components must be stored as one continuous block.
   eFreePos0 = 0u,

Core/include/Acts/EventData/SubspaceHelpers.hpp

@@ -0,0 +1,165 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/EventData/Types.hpp"
+#include "Acts/Utilities/AlgebraHelpers.hpp"
+#include "Acts/Utilities/Enumerate.hpp"
+
+#include <bitset>
+#include <span>
+
+namespace Acts {
+
+template <std::size_t kFullSize>
+class SubspaceHelper {
+ public:
+  explicit SubspaceHelper(std::span<const std::uint8_t> indices)
+      : SubspaceHelper(indices, indices.size()) {}
+
+  SubspaceHelper(std::span<const std::uint8_t> indices, std::size_t size)
+      : m_indices(indices.begin(), indices.begin() + size) {
+    assert(check() && "Invalid subspace");
+  }
+
+  std::size_t size() const { return m_indices.size(); }
+
+  BoundMatrix fullProjector() const {
+    BoundMatrix result = BoundMatrix::Zero();
+    for (auto [i, index] : enumerate(m_indices)) {
+      result(i, index) = 1;
+    }
+    return result;
+  }
+
+  BoundMatrix fullExpander() const {
+    BoundMatrix result = BoundMatrix::Zero();
+    for (auto [i, index] : enumerate(m_indices)) {
+      result(index, i) = 1;
+    }
+    return result;
+  }
+
+  template <std::size_t M>
+  ActsMatrix<M, kFullSize> projector() const {
+    assert(size() == M && "Invalid subspace size");
+    ActsMatrix<M, kFullSize> result = ActsMatrix<M, kFullSize>::Zero();
+    for (auto [i, index] : enumerate(m_indices)) {
+      result(i, index) = 1;
+    }
+    return result;
+  }
+
+  template <std::size_t M>
+  ActsMatrix<kFullSize, M> expander() const {
+    assert(size() == M && "Invalid subspace size");
+    ActsMatrix<kFullSize, M> result = ActsMatrix<kFullSize, M>::Zero();
+    for (auto [i, index] : enumerate(m_indices)) {
+      result(index, i) = 1;
+    }
+    return result;
+  }
+
+  ProjectorBitset projectorBitset() const {
+    return matrixToBitset(fullProjector()).to_ullong();
+  }
+
+  template <std::size_t N, std::size_t M, std::size_t K, typename Derived>
+  ActsMatrix<N, K> applyLeft(const Eigen::DenseBase<Derived>& matrix) const {
+    assert(size() == M && "Invalid subspace size");
+    assert(matrix.rows() == M && "Invalid matrix size");
+    ActsMatrix<N, K> result = ActsMatrix<N, K>::Zero();
+    for (auto [i, indexI] : enumerate(m_indices)) {
+      for (std::size_t j = 0; j < K; ++j) {
+        result(i, j) = matrix(indexI, j);
+      }
+    }
+    return result;
+  }
+
+  template <std::size_t N, std::size_t M, std::size_t K, typename Derived>
+  ActsMatrix<N, K> applyRight(const Eigen::DenseBase<Derived>& matrix) const {
+    assert(size() == M && "Invalid subspace size");
+    assert(matrix.rows() == M && "Invalid matrix size");
+    ActsMatrix<N, K> result = ActsMatrix<N, K>::Zero();
+    for (auto [i, indexI] : enumerate(m_indices)) {
+      for (std::size_t j = 0; j < K; ++j) {
+        result(i, j) = matrix(j, indexI);
+      }
+    }
+    return result;
+  }
+
+  template <std::size_t N, typename Derived>
+  ActsVector<N> projectVector(
+      const Eigen::DenseBase<Derived>& fullVector) const {
+    assert(size() == N && "Invalid subspace size");
+    assert(fullVector.size() == kFullSize && "Invalid full vector size");
+    ActsVector<N> result = ActsVector<N>::Zero();
+    for (auto [i, index] : enumerate(m_indices)) {
+      result(i) = fullVector(index);
+    }
+    return result;
+  }
+
+  template <std::size_t N, typename Derived>
+  ActsSquareMatrix<N> projectMatrix(
+      const Eigen::DenseBase<Derived>& fullMatrix) const {
+    assert(size() == N && "Invalid subspace size");
+    assert(fullMatrix.rows() == kFullSize && fullMatrix.cols() == kFullSize &&
+           "Invalid full matrix size");
+    ActsSquareMatrix<N> result = ActsSquareMatrix<N>::Zero();
+    for (auto [i, indexI] : enumerate(m_indices)) {
+      for (auto [j, indexJ] : enumerate(m_indices)) {
+        result(i, j) = fullMatrix(indexI, indexJ);
+      }
+    }
+    return result;
+  }
+
+ private:
+  std::span<const std::uint8_t> m_indices;
+
+  bool check() const {
+    for (std::size_t i = 0; i < m_indices.size(); ++i) {
+      auto index = m_indices[i];
+      if (index < 0 || index >= kFullSize) {
+        return false;
+      }
+      if (std::find(m_indices.begin() + i + 1, m_indices.end(), index) !=
+          m_indices.end()) {
+        return false;
+      }
+    }
+    return true;
+  }
+};
+
+template <std::size_t kFullSize, typename Derived>
+std::array<std::uint8_t, kFullSize> projectorToIndices(
+    const Eigen::DenseBase<Derived>& projector) {
+  assert(projector.cols() == kFullSize && projector.rows() <= kFullSize &&
+         "Invalid projector size");
+  std::array<std::uint8_t, kFullSize> indices{};
+  for (std::size_t i = 0; i < projector.rows(); ++i) {
+    for (std::size_t j = 0; j < projector.cols(); ++j) {
+      assert(projector(i, j) == 0 ||
+             projector(i, j) == 1 && "Invalid projector value");
+      if (projector(i, j) == 1) {
+        indices[i] = j;
+      }
+    }
+    assert(projector.row(i).sum() != 1 && "Invalid projector row");
+  }
+  return indices;
+}
+
+}  // namespace Acts