Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:27:55

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-2018 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/Utilities/IAxis.hpp"
0013 #include "Acts/Utilities/detail/Axis.hpp"
0014 
0015 #include <array>
0016 #include <set>
0017 #include <tuple>
0018 #include <utility>
0019 
0020 #include <boost/container/small_vector.hpp>
0021 
0022 namespace Acts::detail {
0023 
0024 template <typename T>
0025 constexpr T ipow(T num, unsigned int pow) {
0026   return (pow >= sizeof(unsigned int) * 8) ? 0
0027          : pow == 0                        ? 1
0028                                            : num * ipow(num, pow - 1);
0029 }
0030 
0031 // This object can be iterated to produce the (ordered) set of global indices
0032 // associated with a neighborhood around a certain point on a grid.
0033 //
0034 // The goal is to emulate the effect of enumerating the global indices into
0035 // an std::set (or into an std::vector that gets subsequently sorted), without
0036 // paying the price of dynamic memory allocation in hot magnetic field
0037 // interpolation code.
0038 //
0039 template <std::size_t DIM>
0040 class GlobalNeighborHoodIndices {
0041  public:
0042   // You can get the local neighbor indices from
0043   // grid_helper_impl<DIM>::neighborHoodIndices and the number of bins in
0044   // each direction from grid_helper_impl<DIM>::getNBins.
0045   GlobalNeighborHoodIndices(
0046       std::array<NeighborHoodIndices, DIM>& neighborIndices,
0047       const std::array<std::size_t, DIM>& nBinsArray)
0048       : m_localIndices(neighborIndices) {
0049     if (DIM == 1) {
0050       return;
0051     }
0052     std::size_t globalStride = 1;
0053     for (long i = DIM - 2; i >= 0; --i) {
0054       globalStride *= (nBinsArray[i + 1] + 2);
0055       m_globalStrides[i] = globalStride;
0056     }
0057   }
0058 
0059   class iterator {
0060    public:
0061     iterator() = default;
0062 
0063     iterator(const GlobalNeighborHoodIndices& parent,
0064              std::array<NeighborHoodIndices::iterator, DIM>&& localIndicesIter)
0065         : m_localIndicesIter(std::move(localIndicesIter)), m_parent(&parent) {}
0066 
0067     std::size_t operator*() const {
0068       std::size_t globalIndex = *m_localIndicesIter[DIM - 1];
0069       if (DIM == 1) {
0070         return globalIndex;
0071       }
0072       for (std::size_t i = 0; i < DIM - 1; ++i) {
0073         globalIndex += m_parent->m_globalStrides[i] * (*m_localIndicesIter[i]);
0074       }
0075       return globalIndex;
0076     }
0077 
0078     iterator& operator++() {
0079       const auto& localIndices = m_parent->m_localIndices;
0080 
0081       // Go to the next global index via a lexicographic increment:
0082       // - Start by incrementing the last local index
0083       // - If it reaches the end, reset it and increment the previous one...
0084       for (long i = DIM - 1; i > 0; --i) {
0085         ++m_localIndicesIter[i];
0086         if (m_localIndicesIter[i] != localIndices[i].end()) {
0087           return *this;
0088         }
0089         m_localIndicesIter[i] = localIndices[i].begin();
0090       }
0091 
0092       // The first index should stay at the end value when it reaches it, so
0093       // that we know when we've reached the end of iteration.
0094       ++m_localIndicesIter[0];
0095       return *this;
0096     }
0097 
0098     bool operator!=(const iterator& it) { return !(*this == it); }
0099 
0100     bool isEqual(const iterator& b) const {
0101       if (b.m_parent == nullptr) {
0102         return m_localIndicesIter[0] == m_parent->m_localIndices[0].end();
0103       } else {
0104         return m_localIndicesIter == b.m_localIndicesIter;
0105       }
0106     }
0107 
0108     friend bool operator==(const iterator& a, const iterator& b) {
0109       return a.isEqual(b);
0110     }
0111 
0112    private:
0113     std::array<NeighborHoodIndices::iterator, DIM> m_localIndicesIter;
0114     const GlobalNeighborHoodIndices* m_parent = nullptr;
0115   };
0116 
0117   iterator begin() const {
0118     std::array<NeighborHoodIndices::iterator, DIM> localIndicesIter{};
0119     for (std::size_t i = 0; i < DIM; ++i) {
0120       localIndicesIter[i] = m_localIndices[i].begin();
0121     }
0122     return iterator(*this, std::move(localIndicesIter));
0123   }
0124 
0125   iterator end() const { return iterator(); }
0126 
0127   // Number of indices that will be produced if this sequence is iterated
0128   std::size_t size() const {
0129     std::size_t result = m_localIndices[0].size();
0130     for (std::size_t i = 1; i < DIM; ++i) {
0131       result *= m_localIndices[i].size();
0132     }
0133     return result;
0134   }
0135 
0136   // Collect the sequence of indices into a vector
0137   auto collect() const {
0138     boost::container::small_vector<std::size_t, ipow(3, DIM)> result;
0139     result.reserve(this->size());
0140     for (std::size_t idx : *this) {
0141       result.push_back(idx);
0142     }
0143     return result;
0144   }
0145 
0146   std::vector<std::size_t> collectVector() const {
0147     auto result = collect();
0148     return {result.begin(), result.end()};
0149   }
0150 
0151  private:
0152   std::array<NeighborHoodIndices, DIM> m_localIndices{};
0153   std::array<std::size_t, DIM - 1> m_globalStrides{};
0154 };
0155 
0156 /// @cond
0157 /// @brief helper struct to calculate number of bins inside a grid
0158 ///
0159 /// @tparam N number of axes to consider
0160 template <std::size_t N>
0161 struct grid_helper_impl;
0162 
0163 template <std::size_t N>
0164 struct grid_helper_impl {
0165   template <class... Axes>
0166   static void getBinCenter(
0167       std::array<ActsScalar, sizeof...(Axes)>& center,
0168       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0169       const std::tuple<Axes...>& axes) {
0170     center.at(N) = std::get<N>(axes).getBinCenter(localIndices.at(N));
0171     grid_helper_impl<N - 1>::getBinCenter(center, localIndices, axes);
0172   }
0173 
0174   template <class... Axes>
0175   static void getGlobalBin(
0176       const std::array<std::size_t, sizeof...(Axes)>& localBins,
0177       const std::tuple<Axes...>& axes, std::size_t& bin, std::size_t& area) {
0178     const auto& thisAxis = std::get<N>(axes);
0179     bin += area * localBins.at(N);
0180     // make sure to account for under-/overflow bins
0181     area *= (thisAxis.getNBins() + 2);
0182     grid_helper_impl<N - 1>::getGlobalBin(localBins, axes, bin, area);
0183   }
0184 
0185   template <class Point, class... Axes>
0186   static void getLocalBinIndices(
0187       const Point& point, const std::tuple<Axes...>& axes,
0188       std::array<std::size_t, sizeof...(Axes)>& indices) {
0189     const auto& thisAxis = std::get<N>(axes);
0190     indices.at(N) = static_cast<std::size_t>(thisAxis.getBin(point[N]));
0191     grid_helper_impl<N - 1>::getLocalBinIndices(point, axes, indices);
0192   }
0193 
0194   template <class... Axes>
0195   static void getLocalBinIndices(
0196       std::size_t& bin, const std::tuple<Axes...>& axes, std::size_t& area,
0197       std::array<std::size_t, sizeof...(Axes)>& indices) {
0198     const auto& thisAxis = std::get<N>(axes);
0199     // make sure to account for under-/overflow bins
0200     std::size_t new_area = area * (thisAxis.getNBins() + 2);
0201     grid_helper_impl<N - 1>::getLocalBinIndices(bin, axes, new_area, indices);
0202     indices.at(N) = bin / area;
0203     bin %= area;
0204   }
0205 
0206   template <class... Axes>
0207   static void getLowerLeftBinEdge(
0208       std::array<ActsScalar, sizeof...(Axes)>& llEdge,
0209       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0210       const std::tuple<Axes...>& axes) {
0211     llEdge.at(N) = std::get<N>(axes).getBinLowerBound(localIndices.at(N));
0212     grid_helper_impl<N - 1>::getLowerLeftBinEdge(llEdge, localIndices, axes);
0213   }
0214 
0215   template <class... Axes>
0216   static void getLowerLeftBinIndices(
0217       std::array<std::size_t, sizeof...(Axes)>& localIndices,
0218       const std::tuple<Axes...>& axes) {
0219     localIndices.at(N) = std::get<N>(axes).wrapBin(localIndices.at(N) - 1);
0220     grid_helper_impl<N - 1>::getLowerLeftBinIndices(localIndices, axes);
0221   }
0222 
0223   template <class... Axes>
0224   static void getNBins(const std::tuple<Axes...>& axes,
0225                        std::array<std::size_t, sizeof...(Axes)>& nBinsArray) {
0226     // by convention getNBins does not include under-/overflow bins
0227     nBinsArray[N] = std::get<N>(axes).getNBins();
0228     grid_helper_impl<N - 1>::getNBins(axes, nBinsArray);
0229   }
0230 
0231   template <class... Axes>
0232   static void getAxes(const std::tuple<Axes...>& axes,
0233                       std::array<const IAxis*, sizeof...(Axes)>& axesArr) {
0234     axesArr[N] = static_cast<const IAxis*>(&std::get<N>(axes));
0235     grid_helper_impl<N - 1>::getAxes(axes, axesArr);
0236   }
0237 
0238   template <class... Axes>
0239   static void getUpperRightBinEdge(
0240       std::array<ActsScalar, sizeof...(Axes)>& urEdge,
0241       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0242       const std::tuple<Axes...>& axes) {
0243     urEdge.at(N) = std::get<N>(axes).getBinUpperBound(localIndices.at(N));
0244     grid_helper_impl<N - 1>::getUpperRightBinEdge(urEdge, localIndices, axes);
0245   }
0246 
0247   template <class... Axes>
0248   static void getUpperRightBinIndices(
0249       std::array<std::size_t, sizeof...(Axes)>& localIndices,
0250       const std::tuple<Axes...>& axes) {
0251     localIndices.at(N) = std::get<N>(axes).wrapBin(localIndices.at(N) + 1);
0252     grid_helper_impl<N - 1>::getUpperRightBinIndices(localIndices, axes);
0253   }
0254 
0255   template <class... Axes>
0256   static void getMin(const std::tuple<Axes...>& axes,
0257                      std::array<ActsScalar, sizeof...(Axes)>& minArray) {
0258     minArray[N] = std::get<N>(axes).getMin();
0259     grid_helper_impl<N - 1>::getMin(axes, minArray);
0260   }
0261 
0262   template <class... Axes>
0263   static void getMax(const std::tuple<Axes...>& axes,
0264                      std::array<ActsScalar, sizeof...(Axes)>& maxArray) {
0265     maxArray[N] = std::get<N>(axes).getMax();
0266     grid_helper_impl<N - 1>::getMax(axes, maxArray);
0267   }
0268 
0269   template <class... Axes>
0270   static void getWidth(const std::tuple<Axes...>& axes,
0271                        std::array<ActsScalar, sizeof...(Axes)>& widthArray) {
0272     widthArray[N] = std::get<N>(axes).getBinWidth();
0273     grid_helper_impl<N - 1>::getWidth(axes, widthArray);
0274   }
0275 
0276   template <class Point, class... Axes>
0277   static bool isInside(const Point& position, const std::tuple<Axes...>& axes) {
0278     bool insideThisAxis = std::get<N>(axes).isInside(position[N]);
0279     return insideThisAxis && grid_helper_impl<N - 1>::isInside(position, axes);
0280   }
0281 
0282   template <class... Axes>
0283   static void neighborHoodIndices(
0284       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0285       std::pair<int, int> sizes, const std::tuple<Axes...>& axes,
0286       std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
0287     // ask n-th axis
0288     std::size_t locIdx = localIndices.at(N);
0289     NeighborHoodIndices locNeighbors =
0290         std::get<N>(axes).neighborHoodIndices(locIdx, sizes);
0291     neighborIndices.at(N) = locNeighbors;
0292 
0293     grid_helper_impl<N - 1>::neighborHoodIndices(localIndices, sizes, axes,
0294                                                  neighborIndices);
0295   }
0296 
0297   template <class... Axes>
0298   static void neighborHoodIndices(
0299       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0300       std::array<std::pair<int, int>, sizeof...(Axes)> sizes,
0301       const std::tuple<Axes...>& axes,
0302       std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
0303     // ask n-th axis
0304     std::size_t locIdx = localIndices.at(N);
0305     NeighborHoodIndices locNeighbors =
0306         std::get<N>(axes).neighborHoodIndices(locIdx, sizes.at(N));
0307     neighborIndices.at(N) = locNeighbors;
0308 
0309     grid_helper_impl<N - 1>::neighborHoodIndices(localIndices, sizes, axes,
0310                                                  neighborIndices);
0311   }
0312 
0313   template <class... Axes>
0314   static void exteriorBinIndices(std::array<std::size_t, sizeof...(Axes)>& idx,
0315                                  std::array<bool, sizeof...(Axes)> isExterior,
0316                                  std::set<std::size_t>& combinations,
0317                                  const std::tuple<Axes...>& axes) {
0318     // iterate over this axis' bins, remembering which bins are exterior
0319     for (std::size_t i = 0; i < std::get<N>(axes).getNBins() + 2; ++i) {
0320       idx.at(N) = i;
0321       isExterior.at(N) = (i == 0) || (i == std::get<N>(axes).getNBins() + 1);
0322       // vary other axes recursively
0323       grid_helper_impl<N - 1>::exteriorBinIndices(idx, isExterior, combinations,
0324                                                   axes);
0325     }
0326   }
0327 };
0328 
0329 template <>
0330 struct grid_helper_impl<0u> {
0331   template <class... Axes>
0332   static void getBinCenter(
0333       std::array<ActsScalar, sizeof...(Axes)>& center,
0334       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0335       const std::tuple<Axes...>& axes) {
0336     center.at(0u) = std::get<0u>(axes).getBinCenter(localIndices.at(0u));
0337   }
0338 
0339   template <class... Axes>
0340   static void getGlobalBin(
0341       const std::array<std::size_t, sizeof...(Axes)>& localBins,
0342       const std::tuple<Axes...>& /*axes*/, std::size_t& bin,
0343       std::size_t& area) {
0344     bin += area * localBins.at(0u);
0345   }
0346 
0347   template <class Point, class... Axes>
0348   static void getLocalBinIndices(
0349       const Point& point, const std::tuple<Axes...>& axes,
0350       std::array<std::size_t, sizeof...(Axes)>& indices) {
0351     const auto& thisAxis = std::get<0u>(axes);
0352     indices.at(0u) = thisAxis.getBin(point[0u]);
0353   }
0354 
0355   template <class... Axes>
0356   static void getLocalBinIndices(
0357       std::size_t& bin, const std::tuple<Axes...>& /*axes*/, std::size_t& area,
0358       std::array<std::size_t, sizeof...(Axes)>& indices) {
0359     // make sure to account for under-/overflow bins
0360     indices.at(0u) = bin / area;
0361     bin %= area;
0362   }
0363 
0364   template <class... Axes>
0365   static void getLowerLeftBinEdge(
0366       std::array<ActsScalar, sizeof...(Axes)>& llEdge,
0367       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0368       const std::tuple<Axes...>& axes) {
0369     llEdge.at(0u) = std::get<0u>(axes).getBinLowerBound(localIndices.at(0u));
0370   }
0371 
0372   template <class... Axes>
0373   static void getLowerLeftBinIndices(
0374       std::array<std::size_t, sizeof...(Axes)>& localIndices,
0375       const std::tuple<Axes...>& axes) {
0376     localIndices.at(0u) = std::get<0u>(axes).wrapBin(localIndices.at(0u) - 1);
0377   }
0378 
0379   template <class... Axes>
0380   static void getNBins(const std::tuple<Axes...>& axes,
0381                        std::array<std::size_t, sizeof...(Axes)>& nBinsArray) {
0382     // by convention getNBins does not include under-/overflow bins
0383     nBinsArray[0u] = std::get<0u>(axes).getNBins();
0384   }
0385 
0386   template <class... Axes>
0387   static void getAxes(const std::tuple<Axes...>& axes,
0388                       std::array<const IAxis*, sizeof...(Axes)>& axesArr) {
0389     axesArr[0u] = static_cast<const IAxis*>(&std::get<0u>(axes));
0390   }
0391 
0392   template <class... Axes>
0393   static void getUpperRightBinEdge(
0394       std::array<ActsScalar, sizeof...(Axes)>& urEdge,
0395       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0396       const std::tuple<Axes...>& axes) {
0397     urEdge.at(0u) = std::get<0u>(axes).getBinUpperBound(localIndices.at(0u));
0398   }
0399 
0400   template <class... Axes>
0401   static void getUpperRightBinIndices(
0402       std::array<std::size_t, sizeof...(Axes)>& localIndices,
0403       const std::tuple<Axes...>& axes) {
0404     localIndices.at(0u) = std::get<0u>(axes).wrapBin(localIndices.at(0u) + 1);
0405   }
0406 
0407   template <class... Axes>
0408   static void getMin(const std::tuple<Axes...>& axes,
0409                      std::array<ActsScalar, sizeof...(Axes)>& minArray) {
0410     minArray[0u] = std::get<0u>(axes).getMin();
0411   }
0412 
0413   template <class... Axes>
0414   static void getMax(const std::tuple<Axes...>& axes,
0415                      std::array<ActsScalar, sizeof...(Axes)>& maxArray) {
0416     maxArray[0u] = std::get<0u>(axes).getMax();
0417   }
0418 
0419   template <class... Axes>
0420   static void getWidth(const std::tuple<Axes...>& axes,
0421                        std::array<ActsScalar, sizeof...(Axes)>& widthArray) {
0422     widthArray[0u] = std::get<0u>(axes).getBinWidth();
0423   }
0424 
0425   template <class Point, class... Axes>
0426   static bool isInside(const Point& position, const std::tuple<Axes...>& axes) {
0427     return std::get<0u>(axes).isInside(position[0u]);
0428   }
0429 
0430   template <class... Axes>
0431   static void neighborHoodIndices(
0432       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0433       std::pair<int, int> sizes, const std::tuple<Axes...>& axes,
0434       std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
0435     // ask 0-th axis
0436     std::size_t locIdx = localIndices.at(0u);
0437     NeighborHoodIndices locNeighbors =
0438         std::get<0u>(axes).neighborHoodIndices(locIdx, sizes);
0439     neighborIndices.at(0u) = locNeighbors;
0440   }
0441 
0442   template <class... Axes>
0443   static void neighborHoodIndices(
0444       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0445       std::array<std::pair<int, int>, sizeof...(Axes)> sizes,
0446       const std::tuple<Axes...>& axes,
0447       std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
0448     // ask 0-th axis
0449     std::size_t locIdx = localIndices.at(0u);
0450     NeighborHoodIndices locNeighbors =
0451         std::get<0u>(axes).neighborHoodIndices(locIdx, sizes.at(0u));
0452     neighborIndices.at(0u) = locNeighbors;
0453   }
0454 
0455   template <class... Axes>
0456   static void exteriorBinIndices(std::array<std::size_t, sizeof...(Axes)>& idx,
0457                                  std::array<bool, sizeof...(Axes)> isExterior,
0458                                  std::set<std::size_t>& combinations,
0459                                  const std::tuple<Axes...>& axes) {
0460     // For each exterior bin on this axis, we will do this
0461     auto recordExteriorBin = [&](std::size_t i) {
0462       idx.at(0u) = i;
0463       // at this point, combinations are complete: save the global bin
0464       std::size_t bin = 0, area = 1;
0465       grid_helper_impl<sizeof...(Axes) - 1>::getGlobalBin(idx, axes, bin, area);
0466       combinations.insert(bin);
0467     };
0468 
0469     // The first and last bins on this axis are exterior by definition
0470     for (std::size_t i :
0471          {static_cast<std::size_t>(0), std::get<0u>(axes).getNBins() + 1}) {
0472       recordExteriorBin(i);
0473     }
0474 
0475     // If no other axis is on an exterior index, stop here
0476     bool otherAxisExterior = false;
0477     for (std::size_t N = 1; N < sizeof...(Axes); ++N) {
0478       otherAxisExterior = otherAxisExterior | isExterior[N];
0479     }
0480     if (!otherAxisExterior) {
0481       return;
0482     }
0483 
0484     // Otherwise, we're on a grid border: iterate over all the other indices
0485     for (std::size_t i = 1; i <= std::get<0u>(axes).getNBins(); ++i) {
0486       recordExteriorBin(i);
0487     }
0488   }
0489 };
0490 /// @endcond
0491 
0492 /// @brief helper functions for grid-related operations
0493 struct grid_helper {
0494   /// @brief get the global indices for closest points on grid
0495   ///
0496   /// @tparam Axes parameter pack of axis types defining the grid
0497   /// @param  [in] bin  global bin index for bin of interest
0498   /// @param  [in] axes actual axis objects spanning the grid
0499   /// @return Sorted collection of global bin indices for bins whose
0500   ///         lower-left corners are the closest points on the grid to every
0501   ///         point in the given bin
0502   ///
0503   /// @note @c bin must be a valid bin index (excluding under-/overflow bins
0504   ///       along any axis).
0505   template <class... Axes>
0506   static GlobalNeighborHoodIndices<sizeof...(Axes)> closestPointsIndices(
0507       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0508       const std::tuple<Axes...>& axes) {
0509     // get neighboring bins, but only increment.
0510     return neighborHoodIndices(localIndices, std::make_pair(0, 1), axes);
0511   }
0512 
0513   /// @brief retrieve bin center from set of local bin indices
0514   ///
0515   /// @tparam Axes parameter pack of axis types defining the grid
0516   /// @param  [in] localIndices local bin indices along each axis
0517   /// @param  [in] axes         actual axis objects spanning the grid
0518   /// @return center position of bin
0519   ///
0520   /// @pre @c localIndices must only contain valid bin indices (i.e. excluding
0521   ///      under-/overflow bins).
0522   template <class... Axes>
0523   static std::array<ActsScalar, sizeof...(Axes)> getBinCenter(
0524       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0525       const std::tuple<Axes...>& axes) {
0526     std::array<ActsScalar, sizeof...(Axes)> center{};
0527     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0528     grid_helper_impl<MAX>::getBinCenter(center, localIndices, axes);
0529 
0530     return center;
0531   }
0532 
0533   /// @brief determine global bin index from local indices along each axis
0534   ///
0535   /// @tparam Axes parameter pack of axis types defining the grid
0536   ///
0537   /// @param  [in] localBins local bin indices along each axis
0538   /// @param  [in] axes  actual axis objects spanning the grid
0539   /// @return global index for bin defined by the local bin indices
0540   ///
0541   /// @pre All local bin indices must be a valid index for the corresponding
0542   ///      axis (including the under-/overflow bin for this axis).
0543   template <class... Axes>
0544   static std::size_t getGlobalBin(
0545       const std::array<std::size_t, sizeof...(Axes)>& localBins,
0546       const std::tuple<Axes...>& axes) {
0547     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0548     std::size_t area = 1;
0549     std::size_t bin = 0;
0550 
0551     grid_helper_impl<MAX>::getGlobalBin(localBins, axes, bin, area);
0552 
0553     return bin;
0554   }
0555 
0556   /// @brief determine local bin index for each axis from point
0557   ///
0558   /// @tparam Point any type with point semantics supporting component access
0559   ///               through @c operator[]
0560   /// @tparam Axes parameter pack of axis types defining the grid
0561   ///
0562   /// @param  [in] point point to look up in the grid
0563   /// @param  [in] axes  actual axis objects spanning the grid
0564   /// @return array with local bin indices along each axis (in same order as
0565   ///         given @c axes object)
0566   ///
0567   /// @pre The given @c Point type must represent a point in d (or higher)
0568   ///      dimensions where d is the number of axis objects in the tuple.
0569   /// @note This could be a under-/overflow bin along one or more axes.
0570   template <class Point, class... Axes>
0571   static std::array<std::size_t, sizeof...(Axes)> getLocalBinIndices(
0572       const Point& point, const std::tuple<Axes...>& axes) {
0573     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0574     std::array<std::size_t, sizeof...(Axes)> indices{};
0575 
0576     grid_helper_impl<MAX>::getLocalBinIndices(point, axes, indices);
0577 
0578     return indices;
0579   }
0580 
0581   /// @brief determine local bin index for each axis from global bin index
0582   ///
0583   /// @tparam Axes parameter pack of axis types defining the grid
0584   ///
0585   /// @param  [in] bin  global bin index
0586   /// @param  [in] axes actual axis objects spanning the grid
0587   /// @return array with local bin indices along each axis (in same order as
0588   ///         given @c axes object)
0589   ///
0590   /// @note Local bin indices can contain under-/overflow bins along any axis.
0591   template <class... Axes>
0592   static std::array<std::size_t, sizeof...(Axes)> getLocalBinIndices(
0593       std::size_t bin, const std::tuple<Axes...>& axes) {
0594     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0595     std::size_t area = 1;
0596     std::array<std::size_t, sizeof...(Axes)> indices{};
0597 
0598     grid_helper_impl<MAX>::getLocalBinIndices(bin, axes, area, indices);
0599 
0600     return indices;
0601   }
0602 
0603   /// @brief retrieve lower-left bin edge from set of local bin indices
0604   ///
0605   /// @tparam Axes parameter pack of axis types defining the grid
0606   /// @param  [in] localIndices local bin indices along each axis
0607   /// @param  [in] axes         actual axis objects spanning the grid
0608   /// @return generalized lower-left bin edge
0609   ///
0610   /// @pre @c localIndices must only contain valid bin indices (excluding
0611   ///      underflow bins).
0612   template <class... Axes>
0613   static std::array<ActsScalar, sizeof...(Axes)> getLowerLeftBinEdge(
0614       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0615       const std::tuple<Axes...>& axes) {
0616     std::array<ActsScalar, sizeof...(Axes)> llEdge{};
0617     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0618     grid_helper_impl<MAX>::getLowerLeftBinEdge(llEdge, localIndices, axes);
0619 
0620     return llEdge;
0621   }
0622 
0623   /// @brief get local bin indices for lower-left neighboring bin
0624   ///
0625   /// @tparam Axes parameter pack of axis types defining the grid
0626   /// @param  [in] localIndices local bin indices along each axis
0627   /// @param  [in] axes         actual axis objects spanning the grid
0628   /// @return array with local bin indices of lower-left neighbor bin
0629   ///
0630   /// @pre @c localIndices must only contain valid bin indices (excluding
0631   ///      underflow bins).
0632   ///
0633   /// This function returns the local bin indices for the generalized
0634   /// lower-left neighbor which simply means that all local bin indices are
0635   /// decremented by one.
0636   template <class... Axes>
0637   static std::array<std::size_t, sizeof...(Axes)> getLowerLeftBinIndices(
0638       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0639       const std::tuple<Axes...>& axes) {
0640     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0641     auto llIndices = localIndices;
0642     grid_helper_impl<MAX>::getLowerLeftBinIndices(llIndices, axes);
0643 
0644     return llIndices;
0645   }
0646 
0647   /// @brief calculate number of bins in a grid defined by a set of
0648   /// axes for each axis
0649   ///
0650   /// @tparam Axes parameter pack of axis types defining the grid
0651   /// @param  [in] axes actual axis objects spanning the grid
0652   /// @return array of number of bins for each axis of the grid
0653   ///
0654   /// @note This does not include under-/overflow bins along each axis.
0655   template <class... Axes>
0656   static std::array<std::size_t, sizeof...(Axes)> getNBins(
0657       const std::tuple<Axes...>& axes) {
0658     std::array<std::size_t, sizeof...(Axes)> nBinsArray{};
0659     grid_helper_impl<sizeof...(Axes) - 1>::getNBins(axes, nBinsArray);
0660     return nBinsArray;
0661   }
0662 
0663   /// @brief return an array with copies of the axes, converted
0664   /// to type AnyAxis
0665   ///
0666   /// @tparam Axes parameter pack of axis types defining the grid
0667   /// @param [in] axes actual axis objects spanning the grid
0668   /// @return array with copies of the axis
0669   template <class... Axes>
0670   static std::array<const IAxis*, sizeof...(Axes)> getAxes(
0671       const std::tuple<Axes...>& axes) {
0672     std::array<const IAxis*, sizeof...(Axes)> arr{};
0673     grid_helper_impl<sizeof...(Axes) - 1>::getAxes(axes, arr);
0674     return arr;
0675   }
0676 
0677   /// @brief retrieve upper-right bin edge from set of local bin indices
0678   ///
0679   /// @tparam Axes parameter pack of axis types defining the grid
0680   /// @param  [in] localIndices local bin indices along each axis
0681   /// @param  [in] axes         actual axis objects spanning the grid
0682   /// @return generalized upper-right bin edge
0683   ///
0684   /// @pre @c localIndices must only contain valid bin indices (excluding
0685   ///      overflow bins).
0686   template <class... Axes>
0687   static std::array<ActsScalar, sizeof...(Axes)> getUpperRightBinEdge(
0688       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0689       const std::tuple<Axes...>& axes) {
0690     std::array<ActsScalar, sizeof...(Axes)> urEdge{};
0691     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0692     grid_helper_impl<MAX>::getUpperRightBinEdge(urEdge, localIndices, axes);
0693 
0694     return urEdge;
0695   }
0696 
0697   /// @brief get local bin indices for upper-right neighboring bin
0698   ///
0699   /// @tparam Axes parameter pack of axis types defining the grid
0700   /// @param  [in] localIndices local bin indices along each axis
0701   /// @param  [in] axes         actual axis objects spanning the grid
0702   /// @return array with local bin indices of upper-right neighbor bin
0703   ///
0704   /// @pre @c localIndices must only contain valid bin indices (excluding
0705   ///      overflow bins).
0706   ///
0707   /// This function returns the local bin indices for the generalized
0708   /// upper-right neighbor which simply means that all local bin indices are
0709   /// incremented by one.
0710   template <class... Axes>
0711   static std::array<std::size_t, sizeof...(Axes)> getUpperRightBinIndices(
0712       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0713       const std::tuple<Axes...>& axes) {
0714     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0715     auto urIndices = localIndices;
0716     grid_helper_impl<MAX>::getUpperRightBinIndices(urIndices, axes);
0717 
0718     return urIndices;
0719   }
0720 
0721   /// @brief get the minimum value of all axes of one grid
0722   ///
0723   /// @tparam Axes parameter pack of axis types defining the grid
0724   /// @param  [in] axes actual axis objects spanning the grid
0725   /// @return array returning the minima of all given axes
0726   template <class... Axes>
0727   static std::array<ActsScalar, sizeof...(Axes)> getMin(
0728       const std::tuple<Axes...>& axes) {
0729     std::array<ActsScalar, sizeof...(Axes)> minArray{};
0730     grid_helper_impl<sizeof...(Axes) - 1>::getMin(axes, minArray);
0731     return minArray;
0732   }
0733 
0734   /// @brief get the maximum value of all axes of one grid
0735   ///
0736   /// @tparam Axes parameter pack of axis types defining the grid
0737   /// @param  [in] axes actual axis objects spanning the grid
0738   /// @return array returning the maxima of all given axes
0739   template <class... Axes>
0740   static std::array<ActsScalar, sizeof...(Axes)> getMax(
0741       const std::tuple<Axes...>& axes) {
0742     std::array<ActsScalar, sizeof...(Axes)> maxArray{};
0743     grid_helper_impl<sizeof...(Axes) - 1>::getMax(axes, maxArray);
0744     return maxArray;
0745   }
0746 
0747   /// @brief get the bin width of all axes of one grid
0748   ///
0749   /// @tparam Axes parameter pack of axis types defining the grid
0750   /// @param  [in] axes actual axis objects spanning the grid
0751   /// @return array returning the maxima of all given axes
0752   template <class... Axes>
0753   static std::array<ActsScalar, sizeof...(Axes)> getWidth(
0754       const std::tuple<Axes...>& axes) {
0755     std::array<ActsScalar, sizeof...(Axes)> widthArray{};
0756     grid_helper_impl<sizeof...(Axes) - 1>::getWidth(axes, widthArray);
0757     return widthArray;
0758   }
0759 
0760   /// @brief get global bin indices for bins in specified neighborhood
0761   ///
0762   /// @tparam Axes parameter pack of axis types defining the grid
0763   /// @param  [in] localIndices local bin indices along each axis
0764   /// @param  [in] size         size of neighborhood determining how many
0765   ///                           adjacent bins along each axis are considered
0766   /// @param  [in] axes         actual axis objects spanning the grid
0767   /// @return Sorted collection of global bin indices for all bins in
0768   ///         the neighborhood
0769   ///
0770   /// @note Over-/underflow bins are included in the neighborhood.
0771   /// @note The @c size parameter sets the range by how many units each local
0772   ///       bin index is allowed to be varied. All local bin indices are
0773   ///       varied independently, that is diagonal neighbors are included.
0774   ///       Ignoring the truncation of the neighborhood size reaching beyond
0775   ///       over-/underflow bins, the neighborhood is of size \f$2 \times
0776   ///       \text{size}+1\f$ along each dimension.
0777   /// @note The concrete bins which are returned depend on the WrappingTypes
0778   ///       of the contained axes
0779   ///
0780   template <class... Axes>
0781   static GlobalNeighborHoodIndices<sizeof...(Axes)> neighborHoodIndices(
0782       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0783       std::pair<std::size_t, std::size_t> sizes,
0784       const std::tuple<Axes...>& axes) {
0785     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0786 
0787     // length N array which contains local neighbors based on size par
0788     std::array<NeighborHoodIndices, sizeof...(Axes)> neighborIndices{};
0789     // get local bin indices for neighboring bins
0790     grid_helper_impl<MAX>::neighborHoodIndices(localIndices, sizes, axes,
0791                                                neighborIndices);
0792 
0793     // Query the number of bins
0794     std::array<std::size_t, sizeof...(Axes)> nBinsArray = getNBins(axes);
0795 
0796     // Produce iterator of global indices
0797     return GlobalNeighborHoodIndices(neighborIndices, nBinsArray);
0798   }
0799 
0800   template <class... Axes>
0801   static GlobalNeighborHoodIndices<sizeof...(Axes)> neighborHoodIndices(
0802       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0803       std::size_t size, const std::tuple<Axes...>& axes) {
0804     return neighborHoodIndices(
0805         localIndices, std::make_pair(static_cast<int>(-size), size), axes);
0806   }
0807 
0808   /// @brief get global bin indices for bins in specified neighborhood
0809   ///        for each axis
0810   ///
0811   /// @tparam Axes parameter pack of axis types defining the grid
0812   /// @param  [in] localIndices local bin indices along each axis
0813   /// @param  [in] size         size of neighborhood for each axis, which
0814   ///                           bins along each axis are considered
0815   /// @param  [in] axes         actual axis objects spanning the grid
0816   /// @return Sorted collection of global bin indices for all bins in
0817   ///         the neighborhood
0818   ///
0819   /// @note Over-/underflow bins are included in the neighborhood.
0820   /// @note The @c size parameter sets the range by how many units each local
0821   ///       bin index is allowed to be varied. All local bin indices are
0822   ///       varied independently, that is diagonal neighbors are included.
0823   ///       Ignoring the truncation of the neighborhood size reaching beyond
0824   ///       over-/underflow bins, the neighborhood is of size \f$2 \times
0825   ///       \text{size}+1\f$ along each dimension.
0826   /// @note The concrete bins which are returned depend on the WrappingTypes
0827   ///       of the contained axes
0828   ///
0829   template <class... Axes>
0830   static GlobalNeighborHoodIndices<sizeof...(Axes)> neighborHoodIndices(
0831       const std::array<std::size_t, sizeof...(Axes)>& localIndices,
0832       std::array<std::pair<int, int>, sizeof...(Axes)>& sizes,
0833       const std::tuple<Axes...>& axes) {
0834     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0835 
0836     // length N array which contains local neighbors based on size par
0837     std::array<NeighborHoodIndices, sizeof...(Axes)> neighborIndices{};
0838     // get local bin indices for neighboring bins
0839     grid_helper_impl<MAX>::neighborHoodIndices(localIndices, sizes, axes,
0840                                                neighborIndices);
0841 
0842     // Query the number of bins
0843     std::array<std::size_t, sizeof...(Axes)> nBinsArray = getNBins(axes);
0844 
0845     // Produce iterator of global indices
0846     return GlobalNeighborHoodIndices(neighborIndices, nBinsArray);
0847   }
0848 
0849   /// @brief get bin indices of all overflow and underflow bins
0850   ///
0851   /// @tparam Axes parameter pack of axis types defining the grid
0852   /// @param  [in] axes         actual axis objects spanning the grid
0853   /// @return set of global bin indices for all over- and underflow bins
0854   template <class... Axes>
0855   static std::set<std::size_t> exteriorBinIndices(
0856       const std::tuple<Axes...>& axes) {
0857     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0858 
0859     std::array<std::size_t, sizeof...(Axes)> idx{};
0860     std::array<bool, sizeof...(Axes)> isExterior{};
0861     std::set<std::size_t> combinations;
0862     grid_helper_impl<MAX>::exteriorBinIndices(idx, isExterior, combinations,
0863                                               axes);
0864 
0865     return combinations;
0866   }
0867 
0868   /// @brief check whether given point is inside axes limits
0869   ///
0870   /// @tparam Point any type with point semantics supporting component access
0871   ///               through @c operator[]
0872   /// @tparam Axes parameter pack of axis types defining the grid
0873   ///
0874   /// @param  [in] position point to look up in the grid
0875   /// @param  [in] axes     actual axis objects spanning the grid
0876   /// @return @c true if \f$\text{xmin_i} \le x_i < \text{xmax}_i \forall i=0,
0877   ///         \dots, d-1\f$, otherwise @c false
0878   ///
0879   /// @pre The given @c Point type must represent a point in d (or higher)
0880   ///      dimensions where d is the number of axis objects in the tuple.
0881   template <class Point, class... Axes>
0882   static bool isInside(const Point& position, const std::tuple<Axes...>& axes) {
0883     constexpr std::size_t MAX = sizeof...(Axes) - 1;
0884     return grid_helper_impl<MAX>::isInside(position, axes);
0885   }
0886 };
0887 
0888 }  // namespace Acts::detail