ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Axis.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Axis.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 
11 #include <algorithm>
12 #include <cmath>
13 #include <vector>
14 #include "Acts/Utilities/IAxis.hpp"
16 
17 namespace Acts {
18 
19 namespace detail {
20 
21 // This object can be iterated to produce up to two sequences of integer
22 // indices, corresponding to the half-open integer ranges [begin1, end1[ and
23 // [begin2, end2[.
24 //
25 // The goal is to emulate the effect of enumerating a range of neighbor
26 // indices on an axis (which may go out of bounds and wrap around since we
27 // have AxisBoundaryType::Closed), inserting them into an std::vector, and
28 // discarding duplicates, without paying the price of duplicate removal
29 // and dynamic memory allocation in hot magnetic field interpolation code.
30 //
32  public:
33  NeighborHoodIndices() = default;
34 
35  NeighborHoodIndices(size_t begin, size_t end)
36  : m_begin1(begin), m_end1(end), m_begin2(end), m_end2(end) {}
37 
38  NeighborHoodIndices(size_t begin1, size_t end1, size_t begin2, size_t end2)
39  : m_begin1(begin1), m_end1(end1), m_begin2(begin2), m_end2(end2) {}
40 
41  class iterator {
42  public:
43  iterator() = default;
44 
45  // Specialized constructor for end() iterator
46  iterator(size_t current) : m_current(current), m_wrapped(true) {}
47 
48  iterator(size_t begin1, size_t end1, size_t begin2)
49  : m_current(begin1),
50  m_end1(end1),
51  m_begin2(begin2),
52  m_wrapped(begin1 == begin2) {}
53 
54  size_t operator*() const { return m_current; }
55 
57  ++m_current;
58  if (m_current == m_end1) {
60  m_wrapped = true;
61  }
62  return *this;
63  }
64 
65  bool operator==(const iterator& it) const {
66  return (m_current == it.m_current) && (m_wrapped == it.m_wrapped);
67  }
68 
69  bool operator!=(const iterator& it) const { return !(*this == it); }
70 
71  private:
73  bool m_wrapped;
74  };
75 
76  iterator begin() const { return iterator(m_begin1, m_end1, m_begin2); }
77 
78  iterator end() const { return iterator(m_end2); }
79 
80  // Number of indices that will be produced if this sequence is iterated
81  size_t size() const { return (m_end1 - m_begin1) + (m_end2 - m_begin2); }
82 
83  // Collect the sequence of indices into an std::vector
84  std::vector<size_t> collect() const {
85  std::vector<size_t> result;
86  result.reserve(this->size());
87  for (size_t idx : *this) {
88  result.push_back(idx);
89  }
90  return result;
91  }
92 
93  private:
94  size_t m_begin1 = 0, m_end1 = 0, m_begin2 = 0, m_end2 = 0;
95 };
96 
101 template <AxisBoundaryType bdt>
102 class Axis<AxisType::Equidistant, bdt> final : public IAxis {
103  public:
112  Axis(double xmin, double xmax, size_t nBins)
113  : m_min(xmin),
114  m_max(xmax),
115  m_width((xmax - xmin) / nBins),
116  m_bins(nBins) {}
117 
121  bool isEquidistant() const override { return true; }
122 
126  bool isVariable() const override { return false; }
127 
131  AxisBoundaryType getBoundaryType() const override { return bdt; }
132 
140  NeighborHoodIndices neighborHoodIndices(size_t idx, size_t size = 1) const {
141  return neighborHoodIndices(idx, std::make_pair(size, size));
142  }
143 
154  template <AxisBoundaryType T = bdt,
155  std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
156  NeighborHoodIndices neighborHoodIndices(size_t idx,
157  std::pair<size_t, size_t> sizes = {
158  1, 1}) const {
159  constexpr int min = 0;
160  const int max = getNBins() + 1;
161  const int itmin = std::max(min, int(idx - sizes.first));
162  const int itmax = std::min(max, int(idx + sizes.second));
163  return NeighborHoodIndices(itmin, itmax + 1);
164  }
165 
175  template <AxisBoundaryType T = bdt,
176  std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
177  NeighborHoodIndices neighborHoodIndices(size_t idx,
178  std::pair<size_t, size_t> sizes = {
179  1, 1}) const {
180  if (idx <= 0 || idx >= (getNBins() + 1)) {
181  return NeighborHoodIndices();
182  }
183  constexpr int min = 1;
184  const int max = getNBins();
185  const int itmin = std::max(min, int(idx - sizes.first));
186  const int itmax = std::min(max, int(idx + sizes.second));
187  return NeighborHoodIndices(itmin, itmax + 1);
188  }
189 
199  template <AxisBoundaryType T = bdt,
200  std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
201  NeighborHoodIndices neighborHoodIndices(size_t idx,
202  std::pair<size_t, size_t> sizes = {
203  1, 1}) const {
204  // Handle invalid indices
205  if (idx <= 0 || idx >= (getNBins() + 1)) {
206  return NeighborHoodIndices();
207  }
208 
209  // Handle corner case where user requests more neighbours than the number
210  // of bins on the axis. We do not want to double-count bins in that case.
211  sizes.first %= getNBins();
212  sizes.second %= getNBins();
213  if (sizes.first + sizes.second + 1 > getNBins()) {
214  sizes.second -= (sizes.first + sizes.second + 1) - getNBins();
215  }
216 
217  // If the entire index range is not covered, we must wrap the range of
218  // targeted neighbor indices into the range of valid bin indices. This may
219  // split the range of neighbor indices in two parts:
220  //
221  // Before wraparound - [ XXXXX]XXX
222  // After wraparound - [ XXXX XXXX ]
223  //
224  const int itmin = idx - sizes.first;
225  const int itmax = idx + sizes.second;
226  const size_t itfirst = wrapBin(itmin);
227  const size_t itlast = wrapBin(itmax);
228  if (itfirst <= itlast) {
229  return NeighborHoodIndices(itfirst, itlast + 1);
230  } else {
231  return NeighborHoodIndices(itfirst, getNBins() + 1, 1, itlast + 1);
232  }
233  }
234 
241  template <AxisBoundaryType T = bdt,
242  std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
243  size_t wrapBin(int bin) const {
244  return std::max(std::min(bin, static_cast<int>(getNBins()) + 1), 0);
245  }
246 
253  template <AxisBoundaryType T = bdt,
254  std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
255  size_t wrapBin(int bin) const {
256  return std::max(std::min(bin, static_cast<int>(getNBins())), 1);
257  }
258 
265  template <AxisBoundaryType T = bdt,
266  std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
267  size_t wrapBin(int bin) const {
268  const int w = getNBins();
269  return 1 + (w + ((bin - 1) % w)) % w;
270  // return int(bin<1)*w - int(bin>w)*w + bin;
271  }
272 
283  size_t getBin(double x) const {
284  return wrapBin(std::floor((x - getMin()) / getBinWidth()) + 1);
285  }
286 
290  double getBinWidth(size_t /*bin*/ = 0) const { return m_width; }
291 
302  double getBinLowerBound(size_t bin) const {
303  return getMin() + (bin - 1) * getBinWidth();
304  }
305 
316  double getBinUpperBound(size_t bin) const {
317  return getMin() + bin * getBinWidth();
318  }
319 
327  double getBinCenter(size_t bin) const {
328  return getMin() + (bin - 0.5) * getBinWidth();
329  }
330 
334  double getMax() const override { return m_max; }
335 
339  double getMin() const override { return m_min; }
340 
344  size_t getNBins() const override { return m_bins; }
345 
353  bool isInside(double x) const { return (m_min <= x) && (x < m_max); }
354 
357  std::vector<double> getBinEdges() const override {
358  std::vector<double> binEdges;
359  for (size_t i = 1; i <= m_bins; i++) {
360  binEdges.push_back(getBinLowerBound(i));
361  }
362  binEdges.push_back(getBinUpperBound(m_bins));
363  return binEdges;
364  }
365 
366  private:
368  double m_min;
370  double m_max;
372  double m_width;
374  size_t m_bins;
375 };
376 
381 template <AxisBoundaryType bdt>
382 class Axis<AxisType::Variable, bdt> final : public IAxis {
383  public:
393  Axis(std::vector<double> binEdges) : m_binEdges(std::move(binEdges)) {}
394 
398  bool isEquidistant() const override { return false; }
399 
403  bool isVariable() const override { return true; }
404 
408  AxisBoundaryType getBoundaryType() const override { return bdt; }
409 
417  NeighborHoodIndices neighborHoodIndices(size_t idx, size_t size = 1) const {
418  return neighborHoodIndices(idx, std::make_pair(size, size));
419  }
420 
431  template <AxisBoundaryType T = bdt,
432  std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
433  NeighborHoodIndices neighborHoodIndices(size_t idx,
434  std::pair<size_t, size_t> sizes = {
435  1, 1}) const {
436  constexpr int min = 0;
437  const int max = getNBins() + 1;
438  const int itmin = std::max(min, int(idx - sizes.first));
439  const int itmax = std::min(max, int(idx + sizes.second));
440  return NeighborHoodIndices(itmin, itmax + 1);
441  }
442 
452  template <AxisBoundaryType T = bdt,
453  std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
454  NeighborHoodIndices neighborHoodIndices(size_t idx,
455  std::pair<size_t, size_t> sizes = {
456  1, 1}) const {
457  if (idx <= 0 || idx >= (getNBins() + 1)) {
458  return NeighborHoodIndices();
459  }
460  constexpr int min = 1;
461  const int max = getNBins();
462  const int itmin = std::max(min, int(idx - sizes.first));
463  const int itmax = std::min(max, int(idx + sizes.second));
464  return NeighborHoodIndices(itmin, itmax + 1);
465  }
466 
476  template <AxisBoundaryType T = bdt,
477  std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
478  NeighborHoodIndices neighborHoodIndices(size_t idx,
479  std::pair<size_t, size_t> sizes = {
480  1, 1}) const {
481  // Handle invalid indices
482  if (idx <= 0 || idx >= (getNBins() + 1)) {
483  return NeighborHoodIndices();
484  }
485 
486  // Handle corner case where user requests more neighbours than the number
487  // of bins on the axis. We do not want to double-count bins in that case.
488  sizes.first %= getNBins();
489  sizes.second %= getNBins();
490  if (sizes.first + sizes.second + 1 > getNBins()) {
491  sizes.second -= (sizes.first + sizes.second + 1) - getNBins();
492  }
493 
494  // If the entire index range is not covered, we must wrap the range of
495  // targeted neighbor indices into the range of valid bin indices. This may
496  // split the range of neighbor indices in two parts:
497  //
498  // Before wraparound - [ XXXXX]XXX
499  // After wraparound - [ XXXX XXXX ]
500  //
501  const int itmin = idx - sizes.first;
502  const int itmax = idx + sizes.second;
503  const size_t itfirst = wrapBin(itmin);
504  const size_t itlast = wrapBin(itmax);
505  if (itfirst <= itlast) {
506  return NeighborHoodIndices(itfirst, itlast + 1);
507  } else {
508  return NeighborHoodIndices(itfirst, getNBins() + 1, 1, itlast + 1);
509  }
510  }
511 
518  template <AxisBoundaryType T = bdt,
519  std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
520  size_t wrapBin(int bin) const {
521  return std::max(std::min(bin, static_cast<int>(getNBins()) + 1), 0);
522  }
523 
530  template <AxisBoundaryType T = bdt,
531  std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
532  size_t wrapBin(int bin) const {
533  return std::max(std::min(bin, static_cast<int>(getNBins())), 1);
534  }
535 
542  template <AxisBoundaryType T = bdt,
543  std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
544  size_t wrapBin(int bin) const {
545  const int w = getNBins();
546  return 1 + (w + ((bin - 1) % w)) % w;
547  // return int(bin<1)*w - int(bin>w)*w + bin;
548  }
549 
560  size_t getBin(double x) const {
561  const auto it =
562  std::upper_bound(std::begin(m_binEdges), std::end(m_binEdges), x);
563  return wrapBin(std::distance(std::begin(m_binEdges), it));
564  }
565 
573  double getBinWidth(size_t bin) const {
574  return m_binEdges.at(bin) - m_binEdges.at(bin - 1);
575  }
576 
587  double getBinLowerBound(size_t bin) const { return m_binEdges.at(bin - 1); }
588 
599  double getBinUpperBound(size_t bin) const { return m_binEdges.at(bin); }
600 
608  double getBinCenter(size_t bin) const {
609  return 0.5 * (getBinLowerBound(bin) + getBinUpperBound(bin));
610  }
611 
615  double getMax() const override { return m_binEdges.back(); }
616 
620  double getMin() const override { return m_binEdges.front(); }
621 
625  size_t getNBins() const override { return m_binEdges.size() - 1; }
626 
634  bool isInside(double x) const {
635  return (m_binEdges.front() <= x) && (x < m_binEdges.back());
636  }
637 
640  std::vector<double> getBinEdges() const override { return m_binEdges; }
641 
642  private:
644  std::vector<double> m_binEdges;
645 };
646 } // namespace detail
647 
648 } // namespace Acts