Utilities/detail/MultiAxisHelper.hpp

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