Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-12 07:51:44

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