Acts/Utilities/GridIterator.ipp

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/Utilities/GridIterator.hpp"
0012
0013 #include <numeric>
0014 #include <stdexcept>
0015
0016 namespace Acts {
0017
0018 // Global Iterator
0019 template <typename T, class... Axes>
0020 GridGlobalIterator<T, Axes...>::GridGlobalIterator(const Grid<T, Axes...>& grid,
0021                                                    std::size_t idx)
0022     : m_grid(&grid), m_idx(idx) {}
0023
0024 template <typename T, class... Axes>
0025 GridGlobalIterator<T, Axes...>::GridGlobalIterator(
0026     GridGlobalIterator<T, Axes...>&& other) noexcept
0027     : m_grid(std::exchange(other.m_grid.ptr, nullptr)), m_idx(other.m_idx) {}
0028
0029 template <typename T, class... Axes>
0030 GridGlobalIterator<T, Axes...>& GridGlobalIterator<T, Axes...>::operator=(
0031     GridGlobalIterator<T, Axes...>&& other) noexcept {
0032   m_grid.ptr = std::exchange(other.m_grid.ptr, nullptr);
0033   m_idx = other.m_idx;
0034   return *this;
0035 }
0036
0037 template <typename T, class... Axes>
0038 bool GridGlobalIterator<T, Axes...>::operator==(
0039     const GridGlobalIterator<T, Axes...>& other) const {
0040   // This will always return false if we are comparing two iterators from
0041   // different grids.
0042   // As such a loop from itrStart (from grid A) to itrStop (from grid B) will
0043   // never complete since itrStop will not be reachable from itrStart
0044   return (m_grid.ptr == other.m_grid.ptr) && m_idx == other.m_idx;
0045 }
0046
0047 template <typename T, class... Axes>
0048 auto GridGlobalIterator<T, Axes...>::operator<=>(
0049     const GridGlobalIterator<T, Axes...>& other) const {
0050   // This operator only makes sense if the two iterators we are comparing
0051   // are using the same grid
0052   assert(m_grid.ptr == other.m_grid.ptr);
0053   return m_idx <=> other.m_idx;
0054 }
0055
0056 template <typename T, class... Axes>
0057 GridGlobalIterator<T, Axes...>& GridGlobalIterator<T, Axes...>::operator+=(
0058     const std::size_t offset) {
0059   m_idx += offset;
0060   return *this;
0061 }
0062
0063 template <typename T, class... Axes>
0064 GridGlobalIterator<T, Axes...>& GridGlobalIterator<T, Axes...>::operator-=(
0065     const std::size_t offset) {
0066   m_idx -= offset;
0067   return *this;
0068 }
0069
0070 template <typename T, class... Axes>
0071 GridGlobalIterator<T, Axes...> GridGlobalIterator<T, Axes...>::operator+(
0072     const std::size_t offset) const {
0073   return GridGlobalIterator<T, Axes...>(*m_grid, m_idx + offset);
0074 }
0075
0076 template <typename T, class... Axes>
0077 GridGlobalIterator<T, Axes...> GridGlobalIterator<T, Axes...>::operator-(
0078     const std::size_t offset) const {
0079   return GridGlobalIterator<T, Axes...>(*m_grid, m_idx - offset);
0080 }
0081
0082 template <typename T, class... Axes>
0083 typename GridGlobalIterator<T, Axes...>::difference_type
0084 GridGlobalIterator<T, Axes...>::operator-(
0085     const GridGlobalIterator<T, Axes...>& other) const {
0086   assert(m_grid.ptr == other.m_grid.ptr);
0087   assert(other <= *this);
0088   return m_idx - other.m_idx;
0089 }
0090
0091 template <typename T, class... Axes>
0092 const typename GridGlobalIterator<T, Axes...>::value_type&
0093 GridGlobalIterator<T, Axes...>::operator*() const {
0094   return m_grid->at(m_idx);
0095 }
0096
0097 template <typename T, class... Axes>
0098 GridGlobalIterator<T, Axes...>& GridGlobalIterator<T, Axes...>::operator++() {
0099   ++m_idx;
0100   return *this;
0101 }
0102
0103 template <typename T, class... Axes>
0104 GridGlobalIterator<T, Axes...> GridGlobalIterator<T, Axes...>::operator++(int) {
0105   GridGlobalIterator<T, Axes...> output(*m_grid, m_idx++);
0106   return output;
0107 }
0108
0109 template <typename T, class... Axes>
0110 std::size_t GridGlobalIterator<T, Axes...>::globalBinIndex() const {
0111   return m_idx;
0112 }
0113
0114 template <typename T, class... Axes>
0115 std::array<std::size_t, GridGlobalIterator<T, Axes...>::DIM>
0116 GridGlobalIterator<T, Axes...>::localBinsIndices() const {
0117   return m_grid->localBinsFromGlobalBin(m_idx);
0118 }
0119
0120 // Local Iterator
0121 template <typename T, class... Axes>
0122 GridLocalIterator<T, Axes...>::GridLocalIterator(
0123     const Grid<T, Axes...>& grid, const std::array<std::size_t, DIM>& indices)
0124     : m_grid(&grid),
0125       m_numLocalBins(grid.numLocalBins()),
0126       m_currentIndex(indices) {
0127   // Since the user has not defined a custom navigation pattern, we tell the
0128   // iterator we want to iterate on all the local bins in ascending order from
0129   // 1ul to numLocalBin for that specific axis.
0130   for (std::size_t i(0); i < DIM; ++i) {
0131     m_navigationIndex[i].resize(m_numLocalBins[i]);
0132     std::iota(m_navigationIndex[i].begin(), m_navigationIndex[i].end(), 1ul);
0133   }
0134 }
0135
0136 template <typename T, class... Axes>
0137 GridLocalIterator<T, Axes...>::GridLocalIterator(
0138     const Grid<T, Axes...>& grid, const std::array<std::size_t, DIM>& indices,
0139     std::array<std::vector<std::size_t>, DIM> navigation)
0140     : m_grid(&grid),
0141       m_numLocalBins(grid.numLocalBins()),
0142       m_currentIndex(indices),
0143       m_navigationIndex(std::move(navigation)) {
0144   /// We can allow navigation on only a subset of bins.
0145   /// If the number of specified bins in the navigation for one axis is not
0146   /// zero then override the maximum number of navigation bins instead of using
0147   /// the total number of available bins in the axis
0148   for (std::size_t i(0ul); i < DIM; ++i) {
0149     /// We do not allow empty bin sequences
0150     if (m_navigationIndex[i].size() == 0) {
0151       throw std::invalid_argument(
0152           "Invalid navigation sequence in local grid iterator. No bins "
0153           "specified.");
0154     }
0155     /// Too many bins
0156     if (m_navigationIndex[i].size() > m_numLocalBins[i]) {
0157       throw std::invalid_argument(
0158           "Invalid navigation sequence in local grid iterator. Too many bins "
0159           "specified.");
0160     }
0161     m_numLocalBins[i] = m_navigationIndex[i].size();
0162   }
0163 }
0164
0165 template <typename T, class... Axes>
0166 GridLocalIterator<T, Axes...>::GridLocalIterator(
0167     GridLocalIterator<T, Axes...>&& other) noexcept
0168     : m_grid(std::exchange(other.m_grid.ptr, nullptr)),
0169       m_numLocalBins(other.m_numLocalBins),
0170       m_currentIndex(other.m_currentIndex),
0171       m_navigationIndex(std::move(other.m_navigationIndex)) {}
0172
0173 template <typename T, class... Axes>
0174 GridLocalIterator<T, Axes...>& GridLocalIterator<T, Axes...>::operator=(
0175     GridLocalIterator<T, Axes...>&& other) noexcept {
0176   m_grid.ptr = std::exchange(other.m_grid.ptr, nullptr);
0177   m_numLocalBins = other.m_numLocalBins;
0178   m_currentIndex = other.m_currentIndex;
0179   m_navigationIndex = std::move(other.m_navigationIndex);
0180   return *this;
0181 }
0182
0183 template <typename T, class... Axes>
0184 bool GridLocalIterator<T, Axes...>::operator==(
0185     const GridLocalIterator<T, Axes...>& other) const {
0186   // This will always return false if we are comparing two iterators from
0187   // different grids.
0188   // As such a loop from itrStart (from grid A) to itrStop (from grid B) will
0189   // never complete since itrStop will not be reachable from itrStart
0190   if (m_grid.ptr != other.m_grid.ptr) {
0191     return false;
0192   }
0193
0194   for (std::size_t i(0); i < DIM; ++i) {
0195     if (m_currentIndex[i] != other.m_currentIndex[i]) {
0196       return false;
0197     }
0198   }
0199
0200   return true;
0201 }
0202
0203 template <typename T, class... Axes>
0204 const typename GridLocalIterator<T, Axes...>::value_type&
0205 GridLocalIterator<T, Axes...>::operator*() const {
0206   std::array<std::size_t, DIM> localPositionBin{};
0207   for (std::size_t i(0); i < DIM; ++i) {
0208     localPositionBin[i] = m_navigationIndex[i][m_currentIndex[i]];
0209   }
0210   return m_grid->atLocalBins(localPositionBin);
0211 }
0212
0213 template <typename T, class... Axes>
0214 GridLocalIterator<T, Axes...>& GridLocalIterator<T, Axes...>::operator++() {
0215   increment<DIM - 1>();
0216   return *this;
0217 }
0218
0219 template <typename T, class... Axes>
0220 GridLocalIterator<T, Axes...> GridLocalIterator<T, Axes...>::operator++(int) {
0221   GridLocalIterator<T, Axes...> output(*this);
0222   this->operator++();
0223   return output;
0224 }
0225
0226 template <typename T, class... Axes>
0227 template <std::size_t N>
0228 void GridLocalIterator<T, Axes...>::increment() {
0229   // Check if the current local bin can be incremented, or we reached the end
0230   // of bins in the axis
0231   if (++m_currentIndex[N] < m_numLocalBins[N]) {
0232     return;
0233   }
0234   // We have reached the last bin in the axis, we set the position to 0ul and
0235   // try to increment another axis
0236   if constexpr (N != 0) {
0237     m_currentIndex[N] = 0;
0238     increment<N - 1>();
0239   } else {
0240     m_currentIndex = m_numLocalBins;
0241   }
0242 }
0243
0244 template <typename T, class... Axes>
0245 std::size_t GridLocalIterator<T, Axes...>::globalBinIndex() const {
0246   return m_grid->globalBinFromLocalBins(localBinsIndices());
0247 }
0248
0249 template <typename T, class... Axes>
0250 std::array<std::size_t, GridLocalIterator<T, Axes...>::DIM>
0251 GridLocalIterator<T, Axes...>::localBinsIndices() const {
0252   std::array<std::size_t, DIM> output{};
0253   for (std::size_t i(0); i < DIM; ++i) {
0254     output[i] = m_navigationIndex[i][m_currentIndex[i]];
0255   }
0256   return output;
0257 }
0258
0259 }  // namespace Acts