Skip to content

Commit

Permalink
feat!: changed Grid interface of neighborHoodIndices (#1052)
Browse files Browse the repository at this point in the history
Changed `Grid` interface of `neighborHoodIndices` with differing lower and upper neighbors. This now allows to specify the number of neighbors desired per axis, as opposed to having one number of neighbors for all axes. It also inverts the way neighbors are interpreted:

BREAKING CHANGE: The `Grid` interface changes: f you want 1 neighbor on the left and two on the right, you now have to write `{-1,2}` instead of `{1,2}` as before. This is necessary as this also allows to give "negative" neighbors, i.e. not starting from the current bin, but e.g. only 2 bins left of the current bin: `{-2,-1}`.
  • Loading branch information
robertlangenberg committed Nov 2, 2021
1 parent 8ce4d31 commit 70669b9
Show file tree
Hide file tree
Showing 5 changed files with 197 additions and 51 deletions.
84 changes: 45 additions & 39 deletions Core/include/Acts/Utilities/detail/Axis.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ class Axis<AxisType::Equidistant, bdt> final : public IAxis {
/// @param [in] sizes how many neighboring bins (up/down)
/// @return Set of neighboring bin indices (global)
NeighborHoodIndices neighborHoodIndices(size_t idx, size_t size = 1) const {
return neighborHoodIndices(idx, std::make_pair(size, size));
return neighborHoodIndices(idx, std::make_pair(-size, size));
}

/// @brief Get #size bins which neighbor the one given
Expand All @@ -155,12 +155,13 @@ class Axis<AxisType::Equidistant, bdt> final : public IAxis {
template <AxisBoundaryType T = bdt,
std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
NeighborHoodIndices neighborHoodIndices(size_t idx,
std::pair<size_t, size_t> sizes = {
1, 1}) const {
std::pair<int, int> sizes = {
-1, 1}) const {
constexpr int min = 0;
const int max = getNBins() + 1;
const int itmin = std::max(min, int(idx - sizes.first));
const int itmax = std::min(max, int(idx + sizes.second));
const int itmin = std::clamp(static_cast<int>(idx + sizes.first), min, max);
const int itmax =
std::clamp(static_cast<int>(idx + sizes.second), min, max);
return NeighborHoodIndices(itmin, itmax + 1);
}

Expand All @@ -176,21 +177,22 @@ class Axis<AxisType::Equidistant, bdt> final : public IAxis {
template <AxisBoundaryType T = bdt,
std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
NeighborHoodIndices neighborHoodIndices(size_t idx,
std::pair<size_t, size_t> sizes = {
1, 1}) const {
std::pair<int, int> sizes = {
-1, 1}) const {
if (idx <= 0 || idx >= (getNBins() + 1)) {
return NeighborHoodIndices();
}
constexpr int min = 1;
const int max = getNBins();
const int itmin = std::max(min, int(idx - sizes.first));
const int itmax = std::min(max, int(idx + sizes.second));
const int itmin = std::clamp(static_cast<int>(idx) + sizes.first, min, max);
const int itmax =
std::clamp(static_cast<int>(idx) + sizes.second, min, max);
return NeighborHoodIndices(itmin, itmax + 1);
}

/// @brief Get #size bins which neighbor the one given
///
/// This is the version for Closed
/// This is the version for Closed (i.e. Wrapping)
///
/// @param [in] idx requested bin index
/// @param [in] sizes how many neighboring bins (up/down)
Expand All @@ -200,20 +202,22 @@ class Axis<AxisType::Equidistant, bdt> final : public IAxis {
template <AxisBoundaryType T = bdt,
std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
NeighborHoodIndices neighborHoodIndices(size_t idx,
std::pair<size_t, size_t> sizes = {
1, 1}) const {
std::pair<int, int> sizes = {
-1, 1}) const {
// Handle invalid indices
if (idx <= 0 || idx >= (getNBins() + 1)) {
return NeighborHoodIndices();
}

// Handle corner case where user requests more neighbours than the number
// of bins on the axis. We do not want to ActsScalar-count bins in that
// case.
sizes.first %= getNBins();
sizes.second %= getNBins();
if (sizes.first + sizes.second + 1 > getNBins()) {
sizes.second -= (sizes.first + sizes.second + 1) - getNBins();
// of bins on the axis. All bins are returned in this case.

const int max = getNBins();
sizes.first = std::clamp(sizes.first, -max, max);
sizes.second = std::clamp(sizes.second, -max, max);
if (std::abs(sizes.first - sizes.second) >= max) {
sizes.first = 1 - idx;
sizes.second = max - idx;
}

// If the entire index range is not covered, we must wrap the range of
Expand All @@ -223,14 +227,14 @@ class Axis<AxisType::Equidistant, bdt> final : public IAxis {
// Before wraparound - [ XXXXX]XXX
// After wraparound - [ XXXX XXXX ]
//
const int itmin = idx - sizes.first;
const int itmin = idx + sizes.first;
const int itmax = idx + sizes.second;
const size_t itfirst = wrapBin(itmin);
const size_t itlast = wrapBin(itmax);
if (itfirst <= itlast) {
return NeighborHoodIndices(itfirst, itlast + 1);
} else {
return NeighborHoodIndices(itfirst, getNBins() + 1, 1, itlast + 1);
return NeighborHoodIndices(itfirst, max + 1, 1, itlast + 1);
}
}

Expand Down Expand Up @@ -417,7 +421,7 @@ class Axis<AxisType::Variable, bdt> final : public IAxis {
/// @param [in] size how many neighboring bins
/// @return Set of neighboring bin indices (global)
NeighborHoodIndices neighborHoodIndices(size_t idx, size_t size = 1) const {
return neighborHoodIndices(idx, std::make_pair(size, size));
return neighborHoodIndices(idx, std::make_pair(-size, size));
}

/// @brief Get #size bins which neighbor the one given
Expand All @@ -433,12 +437,12 @@ class Axis<AxisType::Variable, bdt> final : public IAxis {
template <AxisBoundaryType T = bdt,
std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
NeighborHoodIndices neighborHoodIndices(size_t idx,
std::pair<size_t, size_t> sizes = {
1, 1}) const {
std::pair<int, int> sizes = {
-1, 1}) const {
constexpr int min = 0;
const int max = getNBins() + 1;
const int itmin = std::max(min, int(idx - sizes.first));
const int itmax = std::min(max, int(idx + sizes.second));
const int itmin = std::max(min, static_cast<int>(idx) + sizes.first);
const int itmax = std::min(max, static_cast<int>(idx) + sizes.second);
return NeighborHoodIndices(itmin, itmax + 1);
}

Expand All @@ -454,15 +458,15 @@ class Axis<AxisType::Variable, bdt> final : public IAxis {
template <AxisBoundaryType T = bdt,
std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
NeighborHoodIndices neighborHoodIndices(size_t idx,
std::pair<size_t, size_t> sizes = {
1, 1}) const {
std::pair<int, int> sizes = {
-1, 1}) const {
if (idx <= 0 || idx >= (getNBins() + 1)) {
return NeighborHoodIndices();
}
constexpr int min = 1;
const int max = getNBins();
const int itmin = std::max(min, int(idx - sizes.first));
const int itmax = std::min(max, int(idx + sizes.second));
const int itmin = std::max(min, static_cast<int>(idx) + sizes.first);
const int itmax = std::min(max, static_cast<int>(idx) + sizes.second);
return NeighborHoodIndices(itmin, itmax + 1);
}

Expand All @@ -478,20 +482,22 @@ class Axis<AxisType::Variable, bdt> final : public IAxis {
template <AxisBoundaryType T = bdt,
std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
NeighborHoodIndices neighborHoodIndices(size_t idx,
std::pair<size_t, size_t> sizes = {
1, 1}) const {
std::pair<int, int> sizes = {
-1, 1}) const {
// Handle invalid indices
if (idx <= 0 || idx >= (getNBins() + 1)) {
return NeighborHoodIndices();
}

// Handle corner case where user requests more neighbours than the number
// of bins on the axis. We do not want to ActsScalar-count bins in that
// case.
sizes.first %= getNBins();
sizes.second %= getNBins();
if (sizes.first + sizes.second + 1 > getNBins()) {
sizes.second -= (sizes.first + sizes.second + 1) - getNBins();
// of bins on the axis. All bins are returned in this case

const int max = getNBins();
sizes.first = std::clamp(sizes.first, -max, max);
sizes.second = std::clamp(sizes.second, -max, max);
if (std::abs(sizes.first - sizes.second) >= max) {
sizes.first = 1 - idx;
sizes.second = max - idx;
}

// If the entire index range is not covered, we must wrap the range of
Expand All @@ -501,14 +507,14 @@ class Axis<AxisType::Variable, bdt> final : public IAxis {
// Before wraparound - [ XXXXX]XXX
// After wraparound - [ XXXX XXXX ]
//
const int itmin = idx - sizes.first;
const int itmin = idx + sizes.first;
const int itmax = idx + sizes.second;
const size_t itfirst = wrapBin(itmin);
const size_t itlast = wrapBin(itmax);
if (itfirst <= itlast) {
return NeighborHoodIndices(itfirst, itlast + 1);
} else {
return NeighborHoodIndices(itfirst, getNBins() + 1, 1, itlast + 1);
return NeighborHoodIndices(itfirst, max + 1, 1, itlast + 1);
}
}

Expand Down
22 changes: 22 additions & 0 deletions Core/include/Acts/Utilities/detail/Grid.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -404,6 +404,28 @@ class Grid final {
return grid_helper::neighborHoodIndices(localBins, size, m_axes);
}

/// @brief get global bin indices for neighborhood
///
/// @param [in] localBins center bin defined by local bin indices along
/// each axis. If size is negative, center bin
/// is not returned.
/// @param [in] sizePerAxis size of neighborhood for each axis, how many
/// adjacent bins along each axis are considered
/// @return set of global bin indices for all bins in neighborhood
///
/// @note Over-/underflow bins are included in the neighborhood.
/// @note The @c size parameter sets the range by how many units each local
/// bin index is allowed to be varied. All local bin indices are
/// varied independently, that is diagonal neighbors are included.
/// Ignoring the truncation of the neighborhood size reaching beyond
/// over-/underflow bins, the neighborhood is of size \f$2 \times
/// \text{size}+1\f$ along each dimension.
detail::GlobalNeighborHoodIndices<DIM> neighborHoodIndices(
const index_t& localBins,
std::array<std::pair<int, int>, DIM>& sizePerAxis) const {
return grid_helper::neighborHoodIndices(localBins, sizePerAxis, m_axes);
}

/// @brief total number of bins
///
/// @return total number of bins in the grid
Expand Down
77 changes: 74 additions & 3 deletions Core/include/Acts/Utilities/detail/grid_helper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -280,7 +280,7 @@ struct grid_helper_impl {
template <class... Axes>
static void neighborHoodIndices(
const std::array<size_t, sizeof...(Axes)>& localIndices,
std::pair<size_t, size_t> sizes, const std::tuple<Axes...>& axes,
std::pair<int, int> sizes, const std::tuple<Axes...>& axes,
std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
// ask n-th axis
size_t locIdx = localIndices.at(N);
Expand All @@ -292,6 +292,22 @@ struct grid_helper_impl {
neighborIndices);
}

template <class... Axes>
static void neighborHoodIndices(
const std::array<size_t, sizeof...(Axes)>& localIndices,
std::array<std::pair<int, int>, sizeof...(Axes)> sizes,
const std::tuple<Axes...>& axes,
std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
// ask n-th axis
size_t locIdx = localIndices.at(N);
NeighborHoodIndices locNeighbors =
std::get<N>(axes).neighborHoodIndices(locIdx, sizes.at(N));
neighborIndices.at(N) = locNeighbors;

grid_helper_impl<N - 1>::neighborHoodIndices(localIndices, sizes, axes,
neighborIndices);
}

template <class... Axes>
static void exteriorBinIndices(std::array<size_t, sizeof...(Axes)>& idx,
std::array<bool, sizeof...(Axes)> isExterior,
Expand Down Expand Up @@ -412,7 +428,7 @@ struct grid_helper_impl<0u> {
template <class... Axes>
static void neighborHoodIndices(
const std::array<size_t, sizeof...(Axes)>& localIndices,
std::pair<size_t, size_t> sizes, const std::tuple<Axes...>& axes,
std::pair<int, int> sizes, const std::tuple<Axes...>& axes,
std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
// ask 0-th axis
size_t locIdx = localIndices.at(0u);
Expand All @@ -421,6 +437,19 @@ struct grid_helper_impl<0u> {
neighborIndices.at(0u) = locNeighbors;
}

template <class... Axes>
static void neighborHoodIndices(
const std::array<size_t, sizeof...(Axes)>& localIndices,
std::array<std::pair<int, int>, sizeof...(Axes)> sizes,
const std::tuple<Axes...>& axes,
std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices) {
// ask 0-th axis
size_t locIdx = localIndices.at(0u);
NeighborHoodIndices locNeighbors =
std::get<0u>(axes).neighborHoodIndices(locIdx, sizes.at(0u));
neighborIndices.at(0u) = locNeighbors;
}

template <class... Axes>
static void exteriorBinIndices(std::array<size_t, sizeof...(Axes)>& idx,
std::array<bool, sizeof...(Axes)> isExterior,
Expand Down Expand Up @@ -769,7 +798,49 @@ struct grid_helper {
static GlobalNeighborHoodIndices<sizeof...(Axes)> neighborHoodIndices(
const std::array<size_t, sizeof...(Axes)>& localIndices, size_t size,
const std::tuple<Axes...>& axes) {
return neighborHoodIndices(localIndices, std::make_pair(size, size), axes);
return neighborHoodIndices(localIndices, std::make_pair(int(-size), size),
axes);
}

/// @brief get global bin indices for bins in specified neighborhood
/// for each axis
///
/// @tparam Axes parameter pack of axis types defining the grid
/// @param [in] localIndices local bin indices along each axis
/// @param [in] size size of neighborhood for each axis, which
/// bins along each axis are considered
/// @param [in] axes actual axis objects spanning the grid
/// @return Sorted collection of global bin indices for all bins in
/// the neighborhood
///
/// @note Over-/underflow bins are included in the neighborhood.
/// @note The @c size parameter sets the range by how many units each local
/// bin index is allowed to be varied. All local bin indices are
/// varied independently, that is diagonal neighbors are included.
/// Ignoring the truncation of the neighborhood size reaching beyond
/// over-/underflow bins, the neighborhood is of size \f$2 \times
/// \text{size}+1\f$ along each dimension.
/// @note The concrete bins which are returned depend on the WrappingTypes
/// of the contained axes
///
template <class... Axes>
static GlobalNeighborHoodIndices<sizeof...(Axes)> neighborHoodIndices(
const std::array<size_t, sizeof...(Axes)>& localIndices,
std::array<std::pair<int, int>, sizeof...(Axes)>& sizes,
const std::tuple<Axes...>& axes) {
constexpr size_t MAX = sizeof...(Axes) - 1;

// length N array which contains local neighbors based on size par
std::array<NeighborHoodIndices, sizeof...(Axes)> neighborIndices;
// get local bin indices for neighboring bins
grid_helper_impl<MAX>::neighborHoodIndices(localIndices, sizes, axes,
neighborIndices);

// Query the number of bins
std::array<size_t, sizeof...(Axes)> nBinsArray = getNBins(axes);

// Produce iterator of global indices
return GlobalNeighborHoodIndices(neighborIndices, nBinsArray);
}

/// @brief get bin indices of all overflow and underflow bins
Expand Down

0 comments on commit 70669b9

Please sign in to comment.