Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
grid_helper.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file grid_helper.hpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017-2018 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #pragma once
10 
12 #include "Acts/Utilities/IAxis.hpp"
14 
15 #include <array>
16 #include <set>
17 #include <tuple>
18 #include <utility>
19 
20 #include <boost/container/small_vector.hpp>
21 
22 namespace Acts {
23 
24 namespace detail {
25 
26 template <typename T>
27 constexpr T ipow(T num, unsigned int pow) {
28  return (pow >= sizeof(unsigned int) * 8) ? 0
29  : pow == 0 ? 1
30  : num * ipow(num, pow - 1);
31 }
32 
33 // This object can be iterated to produce the (ordered) set of global indices
34 // associated with a neighborhood around a certain point on a grid.
35 //
36 // The goal is to emulate the effect of enumerating the global indices into
37 // an std::set (or into an std::vector that gets subsequently sorted), without
38 // paying the price of dynamic memory allocation in hot magnetic field
39 // interpolation code.
40 //
41 template <size_t DIM>
43  public:
44  // You can get the local neighbor indices from
45  // grid_helper_impl<DIM>::neighborHoodIndices and the number of bins in
46  // each direction from grid_helper_impl<DIM>::getNBins.
48  std::array<NeighborHoodIndices, DIM>& neighborIndices,
49  const std::array<size_t, DIM>& nBinsArray)
50  : m_localIndices(neighborIndices) {
51  if (DIM == 1) {
52  return;
53  }
54  size_t globalStride = 1;
55  for (long i = DIM - 2; i >= 0; --i) {
56  globalStride *= (nBinsArray[i + 1] + 2);
57  m_globalStrides[i] = globalStride;
58  }
59  }
60 
61  class iterator {
62  public:
63  iterator() = default;
64 
66  std::array<NeighborHoodIndices::iterator, DIM>&& localIndicesIter)
67  : m_localIndicesIter(std::move(localIndicesIter)), m_parent(&parent) {}
68 
69  size_t operator*() const {
70  size_t globalIndex = *m_localIndicesIter[DIM - 1];
71  if (DIM == 1) {
72  return globalIndex;
73  }
74  for (size_t i = 0; i < DIM - 1; ++i) {
75  globalIndex += m_parent->m_globalStrides[i] * (*m_localIndicesIter[i]);
76  }
77  return globalIndex;
78  }
79 
81  const auto& localIndices = m_parent->m_localIndices;
82 
83  // Go to the next global index via a lexicographic increment:
84  // - Start by incrementing the last local index
85  // - If it reaches the end, reset it and increment the previous one...
86  for (long i = DIM - 1; i > 0; --i) {
88  if (m_localIndicesIter[i] != localIndices[i].end()) {
89  return *this;
90  }
91  m_localIndicesIter[i] = localIndices[i].begin();
92  }
93 
94  // The first index should stay at the end value when it reaches it, so
95  // that we know when we've reached the end of iteration.
96  ++m_localIndicesIter[0];
97  return *this;
98  }
99 
100  bool operator!=(const iterator& it) { return !(*this == it); }
101 
102  bool isEqual(const iterator& b) const {
103  if (b.m_parent == nullptr) {
104  return m_localIndicesIter[0] == m_parent->m_localIndices[0].end();
105  } else {
107  }
108  }
109 
110  friend bool operator==(const iterator& a, const iterator& b) {
111  return a.isEqual(b);
112  }
113 
114  private:
115  std::array<NeighborHoodIndices::iterator, DIM> m_localIndicesIter;
117  };
118 
119  iterator begin() const {
120  std::array<NeighborHoodIndices::iterator, DIM> localIndicesIter{};
121  for (size_t i = 0; i < DIM; ++i) {
122  localIndicesIter[i] = m_localIndices[i].begin();
123  }
124  return iterator(*this, std::move(localIndicesIter));
125  }
126 
127  iterator end() const { return iterator(); }
128 
129  // Number of indices that will be produced if this sequence is iterated
130  size_t size() const {
131  size_t result = m_localIndices[0].size();
132  for (size_t i = 1; i < DIM; ++i) {
133  result *= m_localIndices[i].size();
134  }
135  return result;
136  }
137 
138  // Collect the sequence of indices into a vector
139  auto collect() const {
140  boost::container::small_vector<size_t, ipow(3, DIM)> result;
141  result.reserve(this->size());
142  for (size_t idx : *this) {
143  result.push_back(idx);
144  }
145  return result;
146  }
147 
148  std::vector<size_t> collectVector() const {
149  auto result = collect();
150  return {result.begin(), result.end()};
151  }
152 
153  private:
154  std::array<NeighborHoodIndices, DIM> m_localIndices{};
155  std::array<size_t, DIM - 1> m_globalStrides{};
156 };
157 
162 template <size_t N>
164 
165 template <size_t N>
167  template <class... Axes>
168  static void getBinCenter(
169  std::array<ActsScalar, sizeof...(Axes)>& center,
170  const std::array<size_t, sizeof...(Axes)>& localIndices,
171  const std::tuple<Axes...>& axes) {
172  center.at(N) = std::get<N>(axes).getBinCenter(localIndices.at(N));
174  }
175 
176  template <class... Axes>
177  static void getGlobalBin(const std::array<size_t, sizeof...(Axes)>& localBins,
178  const std::tuple<Axes...>& axes, size_t& bin,
179  size_t& area) {
180  const auto& thisAxis = std::get<N>(axes);
181  bin += area * localBins.at(N);
182  // make sure to account for under-/overflow bins
183  area *= (thisAxis.getNBins() + 2);
184  grid_helper_impl<N - 1>::getGlobalBin(localBins, axes, bin, area);
185  }
186 
187  template <class Point, class... Axes>
188  static void getLocalBinIndices(const Point& point,
189  const std::tuple<Axes...>& axes,
190  std::array<size_t, sizeof...(Axes)>& indices) {
191  const auto& thisAxis = std::get<N>(axes);
192  indices.at(N) = static_cast<size_t>(thisAxis.getBin(point[N]));
193  grid_helper_impl<N - 1>::getLocalBinIndices(point, axes, indices);
194  }
195 
196  template <class... Axes>
197  static void getLocalBinIndices(size_t& bin, const std::tuple<Axes...>& axes,
198  size_t& area,
199  std::array<size_t, sizeof...(Axes)>& indices) {
200  const auto& thisAxis = std::get<N>(axes);
201  // make sure to account for under-/overflow bins
202  size_t new_area = area * (thisAxis.getNBins() + 2);
203  grid_helper_impl<N - 1>::getLocalBinIndices(bin, axes, new_area, indices);
204  indices.at(N) = bin / area;
205  bin %= area;
206  }
207 
208  template <class... Axes>
209  static void getLowerLeftBinEdge(
210  std::array<ActsScalar, sizeof...(Axes)>& llEdge,
211  const std::array<size_t, sizeof...(Axes)>& localIndices,
212  const std::tuple<Axes...>& axes) {
213  llEdge.at(N) = std::get<N>(axes).getBinLowerBound(localIndices.at(N));
215  }
216 
217  template <class... Axes>
219  std::array<size_t, sizeof...(Axes)>& localIndices,
220  const std::tuple<Axes...>& axes) {
221  localIndices.at(N) = std::get<N>(axes).wrapBin(localIndices.at(N) - 1);
223  }
224 
225  template <class... Axes>
226  static void getNBins(const std::tuple<Axes...>& axes,
227  std::array<size_t, sizeof...(Axes)>& nBinsArray) {
228  // by convention getNBins does not include under-/overflow bins
229  nBinsArray[N] = std::get<N>(axes).getNBins();
230  grid_helper_impl<N - 1>::getNBins(axes, nBinsArray);
231  }
232 
233  template <class... Axes>
234  static void getAxes(const std::tuple<Axes...>& axes,
235  std::array<const IAxis*, sizeof...(Axes)>& axesArr) {
236  axesArr[N] = static_cast<const IAxis*>(&std::get<N>(axes));
237  grid_helper_impl<N - 1>::getAxes(axes, axesArr);
238  }
239 
240  template <class... Axes>
241  static void getUpperRightBinEdge(
242  std::array<ActsScalar, sizeof...(Axes)>& urEdge,
243  const std::array<size_t, sizeof...(Axes)>& localIndices,
244  const std::tuple<Axes...>& axes) {
245  urEdge.at(N) = std::get<N>(axes).getBinUpperBound(localIndices.at(N));
247  }
248 
249  template <class... Axes>
251  std::array<size_t, sizeof...(Axes)>& localIndices,
252  const std::tuple<Axes...>& axes) {
253  localIndices.at(N) = std::get<N>(axes).wrapBin(localIndices.at(N) + 1);
255  }
256 
257  template <class... Axes>
258  static void getMin(const std::tuple<Axes...>& axes,
259  std::array<ActsScalar, sizeof...(Axes)>& minArray) {
260  minArray[N] = std::get<N>(axes).getMin();
261  grid_helper_impl<N - 1>::getMin(axes, minArray);
262  }
263 
264  template <class... Axes>
265  static void getMax(const std::tuple<Axes...>& axes,
266  std::array<ActsScalar, sizeof...(Axes)>& maxArray) {
267  maxArray[N] = std::get<N>(axes).getMax();
268  grid_helper_impl<N - 1>::getMax(axes, maxArray);
269  }
270 
271  template <class... Axes>
272  static void getWidth(const std::tuple<Axes...>& axes,
273  std::array<ActsScalar, sizeof...(Axes)>& widthArray) {
274  widthArray[N] = std::get<N>(axes).getBinWidth();
275  grid_helper_impl<N - 1>::getWidth(axes, widthArray);
276  }
277 
278  template <class Point, class... Axes>
279  static bool isInside(const Point& position, const std::tuple<Axes...>& axes) {
280  bool insideThisAxis = std::get<N>(axes).isInside(position[N]);
281  return insideThisAxis && grid_helper_impl<N - 1>::isInside(position, axes);
282  }
283 
284  template <class... Axes>
285  static void neighborHoodIndices(
286  const std::array<size_t, sizeof...(Axes)>& localIndices,
287  std::pair<int, int> sizes, const std::tuple<Axes...>& axes,
288  std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
289  // ask n-th axis
290  size_t locIdx = localIndices.at(N);
291  NeighborHoodIndices locNeighbors =
292  std::get<N>(axes).neighborHoodIndices(locIdx, sizes);
293  neighborIndices.at(N) = locNeighbors;
294 
296  neighborIndices);
297  }
298 
299  template <class... Axes>
300  static void neighborHoodIndices(
301  const std::array<size_t, sizeof...(Axes)>& localIndices,
302  std::array<std::pair<int, int>, sizeof...(Axes)> sizes,
303  const std::tuple<Axes...>& axes,
304  std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
305  // ask n-th axis
306  size_t locIdx = localIndices.at(N);
307  NeighborHoodIndices locNeighbors =
308  std::get<N>(axes).neighborHoodIndices(locIdx, sizes.at(N));
309  neighborIndices.at(N) = locNeighbors;
310 
312  neighborIndices);
313  }
314 
315  template <class... Axes>
316  static void exteriorBinIndices(std::array<size_t, sizeof...(Axes)>& idx,
317  std::array<bool, sizeof...(Axes)> isExterior,
318  std::set<size_t>& combinations,
319  const std::tuple<Axes...>& axes) {
320  // iterate over this axis' bins, remembering which bins are exterior
321  for (size_t i = 0; i < std::get<N>(axes).getNBins() + 2; ++i) {
322  idx.at(N) = i;
323  isExterior.at(N) = (i == 0) || (i == std::get<N>(axes).getNBins() + 1);
324  // vary other axes recursively
325  grid_helper_impl<N - 1>::exteriorBinIndices(idx, isExterior, combinations,
326  axes);
327  }
328  }
329 };
330 
331 template <>
332 struct grid_helper_impl<0u> {
333  template <class... Axes>
334  static void getBinCenter(
335  std::array<ActsScalar, sizeof...(Axes)>& center,
336  const std::array<size_t, sizeof...(Axes)>& localIndices,
337  const std::tuple<Axes...>& axes) {
338  center.at(0u) = std::get<0u>(axes).getBinCenter(localIndices.at(0u));
339  }
340 
341  template <class... Axes>
342  static void getGlobalBin(const std::array<size_t, sizeof...(Axes)>& localBins,
343  const std::tuple<Axes...>& /*axes*/, size_t& bin,
344  size_t& area) {
345  bin += area * localBins.at(0u);
346  }
347 
348  template <class Point, class... Axes>
349  static void getLocalBinIndices(const Point& point,
350  const std::tuple<Axes...>& axes,
351  std::array<size_t, sizeof...(Axes)>& indices) {
352  const auto& thisAxis = std::get<0u>(axes);
353  indices.at(0u) = thisAxis.getBin(point[0u]);
354  }
355 
356  template <class... Axes>
357  static void getLocalBinIndices(size_t& bin,
358  const std::tuple<Axes...>& /*axes*/,
359  size_t& area,
360  std::array<size_t, sizeof...(Axes)>& indices) {
361  // make sure to account for under-/overflow bins
362  indices.at(0u) = bin / area;
363  bin %= area;
364  }
365 
366  template <class... Axes>
367  static void getLowerLeftBinEdge(
368  std::array<ActsScalar, sizeof...(Axes)>& llEdge,
369  const std::array<size_t, sizeof...(Axes)>& localIndices,
370  const std::tuple<Axes...>& axes) {
371  llEdge.at(0u) = std::get<0u>(axes).getBinLowerBound(localIndices.at(0u));
372  }
373 
374  template <class... Axes>
376  std::array<size_t, sizeof...(Axes)>& localIndices,
377  const std::tuple<Axes...>& axes) {
378  localIndices.at(0u) = std::get<0u>(axes).wrapBin(localIndices.at(0u) - 1);
379  }
380 
381  template <class... Axes>
382  static void getNBins(const std::tuple<Axes...>& axes,
383  std::array<size_t, sizeof...(Axes)>& nBinsArray) {
384  // by convention getNBins does not include under-/overflow bins
385  nBinsArray[0u] = std::get<0u>(axes).getNBins();
386  }
387 
388  template <class... Axes>
389  static void getAxes(const std::tuple<Axes...>& axes,
390  std::array<const IAxis*, sizeof...(Axes)>& axesArr) {
391  axesArr[0u] = static_cast<const IAxis*>(&std::get<0u>(axes));
392  }
393 
394  template <class... Axes>
395  static void getUpperRightBinEdge(
396  std::array<ActsScalar, sizeof...(Axes)>& urEdge,
397  const std::array<size_t, sizeof...(Axes)>& localIndices,
398  const std::tuple<Axes...>& axes) {
399  urEdge.at(0u) = std::get<0u>(axes).getBinUpperBound(localIndices.at(0u));
400  }
401 
402  template <class... Axes>
404  std::array<size_t, sizeof...(Axes)>& localIndices,
405  const std::tuple<Axes...>& axes) {
406  localIndices.at(0u) = std::get<0u>(axes).wrapBin(localIndices.at(0u) + 1);
407  }
408 
409  template <class... Axes>
410  static void getMin(const std::tuple<Axes...>& axes,
411  std::array<ActsScalar, sizeof...(Axes)>& minArray) {
412  minArray[0u] = std::get<0u>(axes).getMin();
413  }
414 
415  template <class... Axes>
416  static void getMax(const std::tuple<Axes...>& axes,
417  std::array<ActsScalar, sizeof...(Axes)>& maxArray) {
418  maxArray[0u] = std::get<0u>(axes).getMax();
419  }
420 
421  template <class... Axes>
422  static void getWidth(const std::tuple<Axes...>& axes,
423  std::array<ActsScalar, sizeof...(Axes)>& widthArray) {
424  widthArray[0u] = std::get<0u>(axes).getBinWidth();
425  }
426 
427  template <class Point, class... Axes>
428  static bool isInside(const Point& position, const std::tuple<Axes...>& axes) {
429  return std::get<0u>(axes).isInside(position[0u]);
430  }
431 
432  template <class... Axes>
433  static void neighborHoodIndices(
434  const std::array<size_t, sizeof...(Axes)>& localIndices,
435  std::pair<int, int> sizes, const std::tuple<Axes...>& axes,
436  std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
437  // ask 0-th axis
438  size_t locIdx = localIndices.at(0u);
439  NeighborHoodIndices locNeighbors =
440  std::get<0u>(axes).neighborHoodIndices(locIdx, sizes);
441  neighborIndices.at(0u) = locNeighbors;
442  }
443 
444  template <class... Axes>
445  static void neighborHoodIndices(
446  const std::array<size_t, sizeof...(Axes)>& localIndices,
447  std::array<std::pair<int, int>, sizeof...(Axes)> sizes,
448  const std::tuple<Axes...>& axes,
449  std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
450  // ask 0-th axis
451  size_t locIdx = localIndices.at(0u);
452  NeighborHoodIndices locNeighbors =
453  std::get<0u>(axes).neighborHoodIndices(locIdx, sizes.at(0u));
454  neighborIndices.at(0u) = locNeighbors;
455  }
456 
457  template <class... Axes>
458  static void exteriorBinIndices(std::array<size_t, sizeof...(Axes)>& idx,
459  std::array<bool, sizeof...(Axes)> isExterior,
460  std::set<size_t>& combinations,
461  const std::tuple<Axes...>& axes) {
462  // For each exterior bin on this axis, we will do this
463  auto recordExteriorBin = [&](size_t i) {
464  idx.at(0u) = i;
465  // at this point, combinations are complete: save the global bin
466  size_t bin = 0, area = 1;
467  grid_helper_impl<sizeof...(Axes) - 1>::getGlobalBin(idx, axes, bin, area);
468  combinations.insert(bin);
469  };
470 
471  // The first and last bins on this axis are exterior by definition
472  for (size_t i :
473  {static_cast<size_t>(0), std::get<0u>(axes).getNBins() + 1}) {
474  recordExteriorBin(i);
475  }
476 
477  // If no other axis is on an exterior index, stop here
478  bool otherAxisExterior = false;
479  for (size_t N = 1; N < sizeof...(Axes); ++N) {
480  otherAxisExterior = otherAxisExterior | isExterior[N];
481  }
482  if (!otherAxisExterior) {
483  return;
484  }
485 
486  // Otherwise, we're on a grid border: iterate over all the other indices
487  for (size_t i = 1; i <= std::get<0u>(axes).getNBins(); ++i) {
488  recordExteriorBin(i);
489  }
490  }
491 };
493 
495 struct grid_helper {
507  template <class... Axes>
509  const std::array<size_t, sizeof...(Axes)>& localIndices,
510  const std::tuple<Axes...>& axes) {
511  // get neighboring bins, but only increment.
512  return neighborHoodIndices(localIndices, std::make_pair(0, 1), axes);
513  }
514 
524  template <class... Axes>
525  static std::array<ActsScalar, sizeof...(Axes)> getBinCenter(
526  const std::array<size_t, sizeof...(Axes)>& localIndices,
527  const std::tuple<Axes...>& axes) {
528  std::array<ActsScalar, sizeof...(Axes)> center{};
529  constexpr size_t MAX = sizeof...(Axes) - 1;
531 
532  return center;
533  }
534 
545  template <class... Axes>
546  static size_t getGlobalBin(
547  const std::array<size_t, sizeof...(Axes)>& localBins,
548  const std::tuple<Axes...>& axes) {
549  constexpr size_t MAX = sizeof...(Axes) - 1;
550  size_t area = 1;
551  size_t bin = 0;
552 
553  grid_helper_impl<MAX>::getGlobalBin(localBins, axes, bin, area);
554 
555  return bin;
556  }
557 
572  template <class Point, class... Axes>
573  static std::array<size_t, sizeof...(Axes)> getLocalBinIndices(
574  const Point& point, const std::tuple<Axes...>& axes) {
575  constexpr size_t MAX = sizeof...(Axes) - 1;
576  std::array<size_t, sizeof...(Axes)> indices{};
577 
578  grid_helper_impl<MAX>::getLocalBinIndices(point, axes, indices);
579 
580  return indices;
581  }
582 
593  template <class... Axes>
594  static std::array<size_t, sizeof...(Axes)> getLocalBinIndices(
595  size_t bin, const std::tuple<Axes...>& axes) {
596  constexpr size_t MAX = sizeof...(Axes) - 1;
597  size_t area = 1;
598  std::array<size_t, sizeof...(Axes)> indices{};
599 
600  grid_helper_impl<MAX>::getLocalBinIndices(bin, axes, area, indices);
601 
602  return indices;
603  }
604 
614  template <class... Axes>
615  static std::array<ActsScalar, sizeof...(Axes)> getLowerLeftBinEdge(
616  const std::array<size_t, sizeof...(Axes)>& localIndices,
617  const std::tuple<Axes...>& axes) {
618  std::array<ActsScalar, sizeof...(Axes)> llEdge{};
619  constexpr size_t MAX = sizeof...(Axes) - 1;
621 
622  return llEdge;
623  }
624 
638  template <class... Axes>
639  static std::array<size_t, sizeof...(Axes)> getLowerLeftBinIndices(
640  const std::array<size_t, sizeof...(Axes)>& localIndices,
641  const std::tuple<Axes...>& axes) {
642  constexpr size_t MAX = sizeof...(Axes) - 1;
643  auto llIndices = localIndices;
645 
646  return llIndices;
647  }
648 
657  template <class... Axes>
658  static std::array<size_t, sizeof...(Axes)> getNBins(
659  const std::tuple<Axes...>& axes) {
660  std::array<size_t, sizeof...(Axes)> nBinsArray{};
661  grid_helper_impl<sizeof...(Axes) - 1>::getNBins(axes, nBinsArray);
662  return nBinsArray;
663  }
664 
671  template <class... Axes>
672  static std::array<const IAxis*, sizeof...(Axes)> getAxes(
673  const std::tuple<Axes...>& axes) {
674  std::array<const IAxis*, sizeof...(Axes)> arr{};
675  grid_helper_impl<sizeof...(Axes) - 1>::getAxes(axes, arr);
676  return arr;
677  }
678 
688  template <class... Axes>
689  static std::array<ActsScalar, sizeof...(Axes)> getUpperRightBinEdge(
690  const std::array<size_t, sizeof...(Axes)>& localIndices,
691  const std::tuple<Axes...>& axes) {
692  std::array<ActsScalar, sizeof...(Axes)> urEdge{};
693  constexpr size_t MAX = sizeof...(Axes) - 1;
695 
696  return urEdge;
697  }
698 
712  template <class... Axes>
713  static std::array<size_t, sizeof...(Axes)> getUpperRightBinIndices(
714  const std::array<size_t, sizeof...(Axes)>& localIndices,
715  const std::tuple<Axes...>& axes) {
716  constexpr size_t MAX = sizeof...(Axes) - 1;
717  auto urIndices = localIndices;
719 
720  return urIndices;
721  }
722 
728  template <class... Axes>
729  static std::array<ActsScalar, sizeof...(Axes)> getMin(
730  const std::tuple<Axes...>& axes) {
731  std::array<ActsScalar, sizeof...(Axes)> minArray{};
732  grid_helper_impl<sizeof...(Axes) - 1>::getMin(axes, minArray);
733  return minArray;
734  }
735 
741  template <class... Axes>
742  static std::array<ActsScalar, sizeof...(Axes)> getMax(
743  const std::tuple<Axes...>& axes) {
744  std::array<ActsScalar, sizeof...(Axes)> maxArray{};
745  grid_helper_impl<sizeof...(Axes) - 1>::getMax(axes, maxArray);
746  return maxArray;
747  }
748 
754  template <class... Axes>
755  static std::array<ActsScalar, sizeof...(Axes)> getWidth(
756  const std::tuple<Axes...>& axes) {
757  std::array<ActsScalar, sizeof...(Axes)> widthArray{};
758  grid_helper_impl<sizeof...(Axes) - 1>::getWidth(axes, widthArray);
759  return widthArray;
760  }
761 
782  template <class... Axes>
784  const std::array<size_t, sizeof...(Axes)>& localIndices,
785  std::pair<size_t, size_t> sizes, const std::tuple<Axes...>& axes) {
786  constexpr size_t MAX = sizeof...(Axes) - 1;
787 
788  // length N array which contains local neighbors based on size par
789  std::array<NeighborHoodIndices, sizeof...(Axes)> neighborIndices{};
790  // get local bin indices for neighboring bins
792  neighborIndices);
793 
794  // Query the number of bins
795  std::array<size_t, sizeof...(Axes)> nBinsArray = getNBins(axes);
796 
797  // Produce iterator of global indices
798  return GlobalNeighborHoodIndices(neighborIndices, nBinsArray);
799  }
800 
801  template <class... Axes>
803  const std::array<size_t, sizeof...(Axes)>& localIndices, size_t size,
804  const std::tuple<Axes...>& axes) {
805  return neighborHoodIndices(localIndices, std::make_pair(int(-size), size),
806  axes);
807  }
808 
830  template <class... Axes>
832  const std::array<size_t, sizeof...(Axes)>& localIndices,
833  std::array<std::pair<int, int>, sizeof...(Axes)>& sizes,
834  const std::tuple<Axes...>& axes) {
835  constexpr size_t MAX = sizeof...(Axes) - 1;
836 
837  // length N array which contains local neighbors based on size par
838  std::array<NeighborHoodIndices, sizeof...(Axes)> neighborIndices{};
839  // get local bin indices for neighboring bins
841  neighborIndices);
842 
843  // Query the number of bins
844  std::array<size_t, sizeof...(Axes)> nBinsArray = getNBins(axes);
845 
846  // Produce iterator of global indices
847  return GlobalNeighborHoodIndices(neighborIndices, nBinsArray);
848  }
849 
855  template <class... Axes>
856  static std::set<size_t> exteriorBinIndices(const std::tuple<Axes...>& axes) {
857  constexpr size_t MAX = sizeof...(Axes) - 1;
858 
859  std::array<size_t, sizeof...(Axes)> idx{};
860  std::array<bool, sizeof...(Axes)> isExterior{};
861  std::set<size_t> combinations;
862  grid_helper_impl<MAX>::exteriorBinIndices(idx, isExterior, combinations,
863  axes);
864 
865  return combinations;
866  }
867 
881  template <class Point, class... Axes>
882  static bool isInside(const Point& position, const std::tuple<Axes...>& axes) {
883  constexpr size_t MAX = sizeof...(Axes) - 1;
884  return grid_helper_impl<MAX>::isInside(position, axes);
885  }
886 };
887 
888 } // namespace detail
889 
890 } // namespace Acts