Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:09:28

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-2023 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 http://mozilla.org/MPL/2.0/.
0008 
0009 #pragma once
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Utilities/IAxis.hpp"
0013 #include "Acts/Utilities/Interpolation.hpp"
0014 #include "Acts/Utilities/detail/grid_helper.hpp"
0015 
0016 #include <array>
0017 #include <numeric>
0018 #include <set>
0019 #include <tuple>
0020 #include <type_traits>
0021 #include <vector>
0022 
0023 namespace Acts {
0024 template <typename T, class... Axes>
0025 class GridGlobalIterator;
0026 
0027 template <typename T, class... Axes>
0028 class GridLocalIterator;
0029 }  // namespace Acts
0030 
0031 namespace Acts {
0032 
0033 /// @brief class for describing a regular multi-dimensional grid
0034 ///
0035 /// @tparam T    type of values stored inside the bins of the grid
0036 /// @tparam Axes parameter pack of axis types defining the grid
0037 ///
0038 /// Class describing a multi-dimensional, regular grid which can store objects
0039 /// in its multi-dimensional bins. Bins are hyper-boxes and can be accessed
0040 /// either by global bin index, local bin indices or position.
0041 ///
0042 /// @note @c T must be default-constructible.
0043 template <typename T, class... Axes>
0044 class Grid final {
0045  public:
0046   /// number of dimensions of the grid
0047   static constexpr std::size_t DIM = sizeof...(Axes);
0048 
0049   /// type of values stored
0050   using value_type = T;
0051   /// reference type to values stored
0052   using reference = value_type&;
0053   /// constant reference type to values stored
0054   using const_reference = const value_type&;
0055   /// type for points in d-dimensional grid space
0056   using point_t = std::array<ActsScalar, DIM>;
0057   /// index type using local bin indices along each axis
0058   using index_t = std::array<std::size_t, DIM>;
0059   /// global iterator type
0060   using global_iterator_t = Acts::GridGlobalIterator<T, Axes...>;
0061   /// local iterator type
0062   using local_iterator_t = Acts::GridLocalIterator<T, Axes...>;
0063 
0064   /// @brief default constructor
0065   ///
0066   /// @param [in] axes actual axis objects spanning the grid
0067   Grid(std::tuple<Axes...>& axes) = delete;
0068 
0069   /// @brief Constructor from const axis tuple, this will allow
0070   /// creating a grid with a different value type from a template
0071   /// grid object.
0072   ///
0073   /// @param axes
0074   Grid(const std::tuple<Axes...>& axes) : m_axes(axes) {
0075     m_values.resize(size());
0076   }
0077 
0078   /// @brief Move constructor from axis tuple
0079   /// @param axes
0080   Grid(std::tuple<Axes...>&& axes) : m_axes(std::move(axes)) {
0081     m_values.resize(size());
0082   }
0083 
0084   /// @brief access value stored in bin for a given point
0085   ///
0086   /// @tparam Point any type with point semantics supporting component access
0087   ///               through @c operator[]
0088   /// @param [in] point point used to look up the corresponding bin in the
0089   ///                   grid
0090   /// @return reference to value stored in bin containing the given point
0091   ///
0092   /// @pre The given @c Point type must represent a point in d (or higher)
0093   ///      dimensions where d is dimensionality of the grid.
0094   ///
0095   /// @note The look-up considers under-/overflow bins along each axis.
0096   ///       Therefore, the look-up will never fail.
0097   //
0098   template <class Point>
0099   reference atPosition(const Point& point) {
0100     return m_values.at(globalBinFromPosition(point));
0101   }
0102 
0103   /// @brief access value stored in bin for a given point
0104   ///
0105   /// @tparam Point any type with point semantics supporting component access
0106   ///               through @c operator[]
0107   /// @param [in] point point used to look up the corresponding bin in the
0108   ///                   grid
0109   /// @return const-reference to value stored in bin containing the given
0110   ///         point
0111   ///
0112   /// @pre The given @c Point type must represent a point in d (or higher)
0113   ///      dimensions where d is dimensionality of the grid.
0114   ///
0115   /// @note The look-up considers under-/overflow bins along each axis.
0116   ///       Therefore, the look-up will never fail.
0117   template <class Point>
0118   const_reference atPosition(const Point& point) const {
0119     return m_values.at(globalBinFromPosition(point));
0120   }
0121 
0122   /// @brief access value stored in bin with given global bin number
0123   ///
0124   /// @param  [in] bin global bin number
0125   /// @return reference to value stored in bin containing the given
0126   ///         point
0127   reference at(std::size_t bin) { return m_values.at(bin); }
0128 
0129   /// @brief access value stored in bin with given global bin number
0130   ///
0131   /// @param  [in] bin global bin number
0132   /// @return const-reference to value stored in bin containing the given
0133   ///         point
0134   const_reference at(std::size_t bin) const { return m_values.at(bin); }
0135 
0136   /// @brief access value stored in bin with given local bin numbers
0137   ///
0138   /// @param  [in] localBins local bin indices along each axis
0139   /// @return reference to value stored in bin containing the given
0140   ///         point
0141   ///
0142   /// @pre All local bin indices must be a valid index for the corresponding
0143   ///      axis (including the under-/overflow bin for this axis).
0144   reference atLocalBins(const index_t& localBins) {
0145     return m_values.at(globalBinFromLocalBins(localBins));
0146   }
0147 
0148   /// @brief access value stored in bin with given local bin numbers
0149   ///
0150   /// @param  [in] localBins local bin indices along each axis
0151   /// @return const-reference to value stored in bin containing the given
0152   ///         point
0153   ///
0154   /// @pre All local bin indices must be a valid index for the corresponding
0155   ///      axis (including the under-/overflow bin for this axis).
0156   const_reference atLocalBins(const index_t& localBins) const {
0157     return m_values.at(globalBinFromLocalBins(localBins));
0158   }
0159 
0160   /// @brief get global bin indices for closest points on grid
0161   ///
0162   /// @tparam Point any type with point semantics supporting component access
0163   ///               through @c operator[]
0164   /// @param [in] position point of interest
0165   /// @return Iterable thatemits the indices of bins whose lower-left corners
0166   ///         are the closest points on the grid to the input.
0167   ///
0168   /// @pre The given @c Point type must represent a point in d (or higher)
0169   ///      dimensions where d is dimensionality of the grid. It must lie
0170   ///      within the grid range (i.e. not within a under-/overflow bin).
0171   template <class Point>
0172   detail::GlobalNeighborHoodIndices<DIM> closestPointsIndices(
0173       const Point& position) const {
0174     return rawClosestPointsIndices(localBinsFromPosition(position));
0175   }
0176 
0177   /// @brief dimensionality of grid
0178   ///
0179   /// @return number of axes spanning the grid
0180   static constexpr std::size_t dimensions() { return DIM; }
0181 
0182   /// @brief get center position of bin with given local bin numbers
0183   ///
0184   /// @param  [in] localBins local bin indices along each axis
0185   /// @return center position of bin
0186   ///
0187   /// @pre All local bin indices must be a valid index for the corresponding
0188   ///      axis (excluding the under-/overflow bins for each axis).
0189   std::array<ActsScalar, DIM> binCenter(const index_t& localBins) const {
0190     return detail::grid_helper::getBinCenter(localBins, m_axes);
0191   }
0192 
0193   /// @brief determine global index for bin containing the given point
0194   ///
0195   /// @tparam Point any type with point semantics supporting component access
0196   ///               through @c operator[]
0197   ///
0198   /// @param  [in] point point to look up in the grid
0199   /// @return global index for bin containing the given point
0200   ///
0201   /// @pre The given @c Point type must represent a point in d (or higher)
0202   ///      dimensions where d is dimensionality of the grid.
0203   /// @note This could be a under-/overflow bin along one or more axes.
0204   template <class Point>
0205   std::size_t globalBinFromPosition(const Point& point) const {
0206     return globalBinFromLocalBins(localBinsFromPosition(point));
0207   }
0208 
0209   /// @brief determine global bin index from local bin indices along each axis
0210   ///
0211   /// @param  [in] localBins local bin indices along each axis
0212   /// @return global index for bin defined by the local bin indices
0213   ///
0214   /// @pre All local bin indices must be a valid index for the corresponding
0215   ///      axis (including the under-/overflow bin for this axis).
0216   std::size_t globalBinFromLocalBins(const index_t& localBins) const {
0217     return detail::grid_helper::getGlobalBin(localBins, m_axes);
0218   }
0219 
0220   /// @brief  determine global bin index of the bin with the lower left edge
0221   ///         closest to the given point for each axis
0222   ///
0223   /// @tparam Point any type with point semantics supporting component access
0224   ///               through @c operator[]
0225   ///
0226   /// @param  [in] point point to look up in the grid
0227   /// @return global index for bin containing the given point
0228   ///
0229   /// @pre The given @c Point type must represent a point in d (or higher)
0230   ///      dimensions where d is dimensionality of the grid.
0231   /// @note This could be a under-/overflow bin along one or more axes.
0232   template <class Point>
0233   std::size_t globalBinFromFromLowerLeftEdge(const Point& point) const {
0234     return globalBinFromLocalBins(localBinsFromLowerLeftEdge(point));
0235   }
0236 
0237   /// @brief  determine local bin index for each axis from the given point
0238   ///
0239   /// @tparam Point any type with point semantics supporting component access
0240   ///               through @c operator[]
0241   ///
0242   /// @param  [in] point point to look up in the grid
0243   /// @return array with local bin indices along each axis (in same order as
0244   ///         given @c axes object)
0245   ///
0246   /// @pre The given @c Point type must represent a point in d (or higher)
0247   ///      dimensions where d is dimensionality of the grid.
0248   /// @note This could be a under-/overflow bin along one or more axes.
0249   template <class Point>
0250   index_t localBinsFromPosition(const Point& point) const {
0251     return detail::grid_helper::getLocalBinIndices(point, m_axes);
0252   }
0253 
0254   /// @brief determine local bin index for each axis from global bin index
0255   ///
0256   /// @param  [in] bin global bin index
0257   /// @return array with local bin indices along each axis (in same order as
0258   ///         given @c axes object)
0259   ///
0260   /// @note Local bin indices can contain under-/overflow bins along the
0261   ///       corresponding axis.
0262   index_t localBinsFromGlobalBin(std::size_t bin) const {
0263     return detail::grid_helper::getLocalBinIndices(bin, m_axes);
0264   }
0265 
0266   /// @brief  determine local bin index of the bin with the lower left edge
0267   ///         closest to the given point for each axis
0268   ///
0269   /// @tparam Point any type with point semantics supporting component access
0270   ///               through @c operator[]
0271   ///
0272   /// @param  [in] point point to look up in the grid
0273   /// @return array with local bin indices along each axis (in same order as
0274   ///         given @c axes object)
0275   ///
0276   /// @pre The given @c Point type must represent a point in d (or higher)
0277   ///      dimensions where d is dimensionality of the grid.
0278   /// @note This could be a under-/overflow bin along one or more axes.
0279   template <class Point>
0280   index_t localBinsFromLowerLeftEdge(const Point& point) const {
0281     Point shiftedPoint;
0282     point_t width = detail::grid_helper::getWidth(m_axes);
0283     for (std::size_t i = 0; i < DIM; i++) {
0284       shiftedPoint[i] = point[i] + width[i] / 2;
0285     }
0286     return detail::grid_helper::getLocalBinIndices(shiftedPoint, m_axes);
0287   }
0288 
0289   /// @brief retrieve lower-left bin edge from set of local bin indices
0290   ///
0291   /// @param  [in] localBins local bin indices along each axis
0292   /// @return generalized lower-left bin edge position
0293   ///
0294   /// @pre @c localBins must only contain valid bin indices (excluding
0295   ///      underflow bins).
0296   point_t lowerLeftBinEdge(const index_t& localBins) const {
0297     return detail::grid_helper::getLowerLeftBinEdge(localBins, m_axes);
0298   }
0299 
0300   /// @brief retrieve upper-right bin edge from set of local bin indices
0301   ///
0302   /// @param  [in] localBins local bin indices along each axis
0303   /// @return generalized upper-right bin edge position
0304   ///
0305   /// @pre @c localBins must only contain valid bin indices (excluding
0306   ///      overflow bins).
0307   point_t upperRightBinEdge(const index_t& localBins) const {
0308     return detail::grid_helper::getUpperRightBinEdge(localBins, m_axes);
0309   }
0310 
0311   /// @brief get bin width along each specific axis
0312   ///
0313   /// @return array giving the bin width alonf all axes
0314   point_t binWidth() const { return detail::grid_helper::getWidth(m_axes); }
0315 
0316   /// @brief get number of bins along each specific axis
0317   ///
0318   /// @return array giving the number of bins along all axes
0319   ///
0320   /// @note Not including under- and overflow bins
0321   index_t numLocalBins() const { return detail::grid_helper::getNBins(m_axes); }
0322 
0323   /// @brief get the minimum value of all axes of one grid
0324   ///
0325   /// @return array returning the minima of all given axes
0326   point_t minPosition() const { return detail::grid_helper::getMin(m_axes); }
0327 
0328   /// @brief get the maximum value of all axes of one grid
0329   ///
0330   /// @return array returning the maxima of all given axes
0331   point_t maxPosition() const { return detail::grid_helper::getMax(m_axes); }
0332 
0333   /// @brief set all overflow and underflow bins to a certain value
0334   ///
0335   /// @param [in] value value to be inserted in every overflow and underflow
0336   ///                   bin of the grid.
0337   ///
0338   void setExteriorBins(const value_type& value) {
0339     for (std::size_t index : detail::grid_helper::exteriorBinIndices(m_axes)) {
0340       at(index) = value;
0341     }
0342   }
0343 
0344   /// @brief interpolate grid values to given position
0345   ///
0346   /// @tparam Point type specifying geometric positions
0347   /// @tparam U     dummy template parameter identical to @c T
0348   ///
0349   /// @param [in] point location to which to interpolate grid values. The
0350   ///                   position must be within the grid dimensions and not
0351   ///                   lie in an under-/overflow bin along any axis.
0352   ///
0353   /// @return interpolated value at given position
0354   ///
0355   /// @pre The given @c Point type must represent a point in d (or higher)
0356   ///      dimensions where d is dimensionality of the grid.
0357   ///
0358   /// @note This function is available only if the following conditions are
0359   /// fulfilled:
0360   /// - Given @c U and @c V of value type @c T as well as two @c ActsScalar
0361   /// @c a and @c b, then the following must be a valid expression <tt>a * U + b
0362   /// * V</tt> yielding an object which is (implicitly) convertible to @c T.
0363   /// - @c Point must represent a d-dimensional position and support
0364   /// coordinate access using @c operator[] which should return a @c
0365   /// ActsScalar (or a value which is implicitly convertible). Coordinate
0366   /// indices must start at 0.
0367   /// @note Bin values are interpreted as being the field values at the
0368   /// lower-left corner of the corresponding hyper-box.
0369   template <class Point, typename U = T,
0370             typename = std::enable_if_t<
0371                 detail::can_interpolate<Point, std::array<ActsScalar, DIM>,
0372                                         std::array<ActsScalar, DIM>, U>::value>>
0373   T interpolate(const Point& point) const {
0374     // there are 2^DIM corner points used during the interpolation
0375     constexpr std::size_t nCorners = 1 << DIM;
0376 
0377     // construct vector of pairs of adjacent bin centers and values
0378     std::array<value_type, nCorners> neighbors{};
0379 
0380     // get local indices for current bin
0381     // value of bin is interpreted as being the field value at its lower left
0382     // corner
0383     const auto& llIndices = localBinsFromPosition(point);
0384 
0385     // get global indices for all surrounding corner points
0386     const auto& closestIndices = rawClosestPointsIndices(llIndices);
0387 
0388     // get values on grid points
0389     std::size_t i = 0;
0390     for (std::size_t index : closestIndices) {
0391       neighbors.at(i++) = at(index);
0392     }
0393 
0394     return Acts::interpolate(point, lowerLeftBinEdge(llIndices),
0395                              upperRightBinEdge(llIndices), neighbors);
0396   }
0397 
0398   /// @brief check whether given point is inside grid limits
0399   ///
0400   /// @return @c true if \f$\text{xmin_i} \le x_i < \text{xmax}_i \forall i=0,
0401   ///         \dots, d-1\f$, otherwise @c false
0402   ///
0403   /// @pre The given @c Point type must represent a point in d (or higher)
0404   ///      dimensions where d is dimensionality of the grid.
0405   ///
0406   /// @post If @c true is returned, the global bin containing the given point
0407   ///       is a valid bin, i.e. it is neither a underflow nor an overflow bin
0408   ///       along any axis.
0409   template <class Point>
0410   bool isInside(const Point& position) const {
0411     return detail::grid_helper::isInside(position, m_axes);
0412   }
0413 
0414   /// @brief get global bin indices for neighborhood
0415   ///
0416   /// @param [in] localBins center bin defined by local bin indices along each
0417   ///                       axis
0418   /// @param [in] size      size of neighborhood determining how many adjacent
0419   ///                       bins along each axis are considered
0420   /// @return set of global bin indices for all bins in neighborhood
0421   ///
0422   /// @note Over-/underflow bins are included in the neighborhood.
0423   /// @note The @c size parameter sets the range by how many units each local
0424   ///       bin index is allowed to be varied. All local bin indices are
0425   ///       varied independently, that is diagonal neighbors are included.
0426   ///       Ignoring the truncation of the neighborhood size reaching beyond
0427   ///       over-/underflow bins, the neighborhood is of size \f$2 \times
0428   ///       \text{size}+1\f$ along each dimension.
0429   detail::GlobalNeighborHoodIndices<DIM> neighborHoodIndices(
0430       const index_t& localBins, std::size_t size = 1u) const {
0431     return detail::grid_helper::neighborHoodIndices(localBins, size, m_axes);
0432   }
0433 
0434   /// @brief get global bin   indices for neighborhood
0435   ///
0436   /// @param [in] localBins   center bin defined by local bin indices along
0437   ///                         each axis. If size is negative, center bin
0438   ///                         is not returned.
0439   /// @param [in] sizePerAxis size of neighborhood for each axis, how many
0440   ///                         adjacent bins along each axis are considered
0441   /// @return set of global bin indices for all bins in neighborhood
0442   ///
0443   /// @note Over-/underflow bins are included in the neighborhood.
0444   /// @note The @c size parameter sets the range by how many units each local
0445   ///       bin index is allowed to be varied. All local bin indices are
0446   ///       varied independently, that is diagonal neighbors are included.
0447   ///       Ignoring the truncation of the neighborhood size reaching beyond
0448   ///       over-/underflow bins, the neighborhood is of size \f$2 \times
0449   ///       \text{size}+1\f$ along each dimension.
0450   detail::GlobalNeighborHoodIndices<DIM> neighborHoodIndices(
0451       const index_t& localBins,
0452       std::array<std::pair<int, int>, DIM>& sizePerAxis) const {
0453     return detail::grid_helper::neighborHoodIndices(localBins, sizePerAxis,
0454                                                     m_axes);
0455   }
0456 
0457   /// @brief total number of bins
0458   ///
0459   /// @return total number of bins in the grid
0460   ///
0461   /// @note This number contains under-and overflow bins along all axes.
0462   std::size_t size(bool fullCounter = true) const {
0463     index_t nBinsArray = numLocalBins();
0464     std::size_t current_size = 1;
0465     // add under-and overflow bins for each axis and multiply all bins
0466     if (fullCounter) {
0467       for (const auto& value : nBinsArray) {
0468         current_size *= value + 2;
0469       }
0470     }
0471     // ignore under-and overflow bins for each axis and multiply all bins
0472     else {
0473       for (const auto& value : nBinsArray) {
0474         current_size *= value;
0475       }
0476     }
0477     return current_size;
0478   }
0479 
0480   /// @brief Convenience function to convert the type of the grid
0481   /// to hold another object type.
0482   ///
0483   /// @tparam U the new grid value type
0484   ///
0485   /// @return a new grid with the same axes and a different value type
0486   template <typename U>
0487   Grid<U, Axes...> convertType() const {
0488     Grid<U, Axes...> cGrid(m_axes);
0489     return cGrid;
0490   }
0491 
0492   /// @brief Convenience function to convert the type of the grid
0493   /// to hold another object type.
0494   ///
0495   /// @tparam converter_t the converter type
0496   ///
0497   /// This is designed to be most flexible with a converter object
0498   /// as a visitor. If needed, such a visitor could also use
0499   /// caching or other techniques to speed up the conversion.
0500   ///
0501   /// @param cVisitor the converter object as visitor
0502   ///
0503   /// @return a new grid with the same axes and a different value type
0504   template <typename converter_t>
0505   Grid<typename converter_t::value_type, Axes...> convertGrid(
0506       converter_t& cVisitor) const {
0507     Grid<typename converter_t::value_type, Axes...> cGrid(m_axes);
0508     // Loop through the values and convert them
0509     for (std::size_t i = 0; i < size(); i++) {
0510       cGrid.at(i) = cVisitor(at(i));
0511     }
0512     return cGrid;
0513   }
0514 
0515   /// @brief get the axes as a tuple
0516   const std::tuple<Axes...>& axesTuple() const { return m_axes; }
0517 
0518   /// @brief get the axes as an array of IAxis pointers
0519   std::array<const IAxis*, DIM> axes() const {
0520     return detail::grid_helper::getAxes(m_axes);
0521   }
0522 
0523   /// begin iterator for global bins
0524   global_iterator_t begin() const { return global_iterator_t(*this, 0); }
0525 
0526   /// end iterator for global bins
0527   global_iterator_t end() const { return global_iterator_t(*this, size()); }
0528 
0529   /// @brief begin iterator for local bins
0530   ///
0531   /// @param navigator is local navigator for the grid
0532   local_iterator_t begin(
0533       const std::array<std::vector<std::size_t>, DIM>& navigator) const {
0534     std::array<std::size_t, DIM> localBin{};
0535     return local_iterator_t(*this, std::move(localBin), navigator);
0536   }
0537 
0538   /// @brief end iterator for local bins
0539   ///
0540   /// @param navigator is local navigator for the grid
0541   local_iterator_t end(
0542       const std::array<std::vector<std::size_t>, DIM>& navigator) const {
0543     std::array<std::size_t, DIM> endline{};
0544     for (std::size_t i(0ul); i < DIM; ++i) {
0545       endline[i] = navigator[i].size();
0546     }
0547     return local_iterator_t(*this, std::move(endline), navigator);
0548   }
0549 
0550  private:
0551   /// set of axis defining the multi-dimensional grid
0552   std::tuple<Axes...> m_axes;
0553   /// linear value store for each bin
0554   std::vector<T> m_values;
0555 
0556   // Part of closestPointsIndices that goes after local bins resolution.
0557   // Used as an interpolation performance optimization, but not exposed as it
0558   // doesn't make that much sense from an API design standpoint.
0559   detail::GlobalNeighborHoodIndices<DIM> rawClosestPointsIndices(
0560       const index_t& localBins) const {
0561     return detail::grid_helper::closestPointsIndices(localBins, m_axes);
0562   }
0563 };
0564 
0565 }  // namespace Acts