Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:11:08

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