Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-11-07 09:15:24

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