Subspace.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Subspace.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2020 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 <array>
#include <cassert>
#include <cstddef>
#include <cstdint>

namespace Acts {
namespace detail {


template <size_t kFullSize, size_t kSize>
class FixedSizeSubspace {
  static_assert(kFullSize <= static_cast<size_t>(UINT8_MAX),
                "Full vector space size is larger than the supported range");
  static_assert(1u <= kSize, "Subspace size must be at least 1");
  static_assert(kSize <= kFullSize,
                "Subspace can only be as large as the full space");

  template <typename source_t>
  using SubspaceVectorFor = Eigen::Matrix<typename source_t::Scalar, kSize, 1>;
  template <typename source_t>
  using FullspaceVectorFor =
      Eigen::Matrix<typename source_t::Scalar, kFullSize, 1>;
  template <typename scalar_t>
  using ProjectionMatrix = Eigen::Matrix<scalar_t, kSize, kFullSize>;
  template <typename scalar_t>
  using ExpansionMatrix = Eigen::Matrix<scalar_t, kFullSize, kSize>;

  // the functionality could also be implemented using a std::bitset where each
  // bit corresponds to an axis in the fullspace and set bits indicate which
  // bits make up the subspace. this would be a more compact representation but
  // complicates the implementation since we can not easily iterate over the
  // indices of the subspace. storing the subspace indices directly requires a
  // bit more memory but is easier to work with. for our typical use cases with
  // n<=8, this still takes only 64bit of memory.
  std::array<uint8_t, kSize> m_axes;

 public:
  template <typename index_t>
  constexpr FixedSizeSubspace(const std::array<index_t, kSize>& indices) {
    for (size_t i = 0u; i < kSize; ++i) {
      assert((indices[i] < kFullSize) and
             "Axis indices must be within the full space");
      if (0u < i) {
        assert((indices[i - 1u] < indices[i]) and
               "Axis indices must be unique and ordered");
      }
    }
    for (size_t i = 0; i < kSize; ++i) {
      m_axes[i] = static_cast<uint8_t>(indices[i]);
    }
  }
  // The subset can not be constructed w/o defining its axis indices.
  FixedSizeSubspace() = delete;
  FixedSizeSubspace(const FixedSizeSubspace&) = default;
  FixedSizeSubspace(FixedSizeSubspace&&) = default;
  FixedSizeSubspace& operator=(const FixedSizeSubspace&) = default;
  FixedSizeSubspace& operator=(FixedSizeSubspace&&) = default;

  static constexpr size_t size() { return kSize; }
  static constexpr size_t fullSize() { return kFullSize; }

  constexpr const std::array<uint8_t, kSize>& indices() const { return m_axes; }

  constexpr bool contains(size_t index) const {
    bool isContained = false;
    // always iterate over all elements to avoid branching and hope the compiler
    // can optimise this for us.
    for (auto a : m_axes) {
      isContained = (isContained or (a == index));
    }
    return isContained;
  }

  template <typename fullspace_vector_t>
  auto projectVector(const Eigen::MatrixBase<fullspace_vector_t>& full) const
      -> SubspaceVectorFor<fullspace_vector_t> {
    EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(fullspace_vector_t, kFullSize);

    SubspaceVectorFor<fullspace_vector_t> sub;
    for (auto i = 0u; i < kSize; ++i) {
      sub[i] = full[m_axes[i]];
    }
    return sub;
  }

  template <typename subspace_vector_t>
  auto expandVector(const Eigen::MatrixBase<subspace_vector_t>& sub) const
      -> FullspaceVectorFor<subspace_vector_t> {
    EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(subspace_vector_t, kSize);

    FullspaceVectorFor<subspace_vector_t> full;
    full.setZero();
    for (auto i = 0u; i < kSize; ++i) {
      full[m_axes[i]] = sub[i];
    }
    return full;
  }

  template <typename scalar_t>
  auto projector() const -> ProjectionMatrix<scalar_t> {
    ProjectionMatrix<scalar_t> proj;
    proj.setZero();
    for (auto i = 0u; i < kSize; ++i) {
      proj(i, m_axes[i]) = 1;
    }
    return proj;
  }

  template <typename scalar_t>
  auto expander() const -> ExpansionMatrix<scalar_t> {
    ExpansionMatrix<scalar_t> expn;
    expn.setZero();
    for (auto i = 0u; i < kSize; ++i) {
      expn(m_axes[i], i) = 1;
    }
    return expn;
  }
};


}  // namespace detail
}  // namespace Acts