Grid.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Grid.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2017-2018 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/Utilities/detail/GridFwd.hpp"

#include <array>
#include <numeric>
#include <set>
#include <tuple>
#include <type_traits>
#include <vector>
#include "Acts/Utilities/IAxis.hpp"
#include "Acts/Utilities/Interpolation.hpp"
#include "Acts/Utilities/detail/grid_helper.hpp"

namespace Acts {

namespace detail {

template <typename T, class... Axes>
class Grid final {
 public:
  static constexpr size_t DIM = sizeof...(Axes);

  using value_type = T;
  using reference = value_type&;
  using const_reference = const value_type&;
  using point_t = std::array<double, DIM>;
  using index_t = std::array<size_t, DIM>;

  Grid(std::tuple<Axes...> axes) : m_axes(std::move(axes)) {
    m_values.resize(size());
  }

  //
  template <class Point>
  reference atPosition(const Point& point) {
    return m_values.at(globalBinFromPosition(point));
  }

  template <class Point>
  const_reference atPosition(const Point& point) const {
    return m_values.at(globalBinFromPosition(point));
  }

  reference at(size_t bin) { return m_values.at(bin); }

  const_reference at(size_t bin) const { return m_values.at(bin); }

  reference atLocalBins(const index_t& localBins) {
    return m_values.at(globalBinFromLocalBins(localBins));
  }

  const_reference atLocalBins(const index_t& localBins) const {
    return m_values.at(globalBinFromLocalBins(localBins));
  }

  template <class Point>
  detail::GlobalNeighborHoodIndices<DIM> closestPointsIndices(
      const Point& position) const {
    return rawClosestPointsIndices(localBinsFromPosition(position));
  }

  static constexpr size_t dimensions() { return DIM; }

  std::array<double, DIM> binCenter(const index_t& localBins) const {
    return grid_helper::getBinCenter(localBins, m_axes);
  }

  template <class Point>
  size_t globalBinFromPosition(const Point& point) const {
    return globalBinFromLocalBins(localBinsFromPosition(point));
  }

  size_t globalBinFromLocalBins(const index_t& localBins) const {
    return grid_helper::getGlobalBin(localBins, m_axes);
  }

  template <class Point>
  size_t globalBinFromFromLowerLeftEdge(const Point& point) const {
    return globalBinFromLocalBins(localBinsFromLowerLeftEdge(point));
  }

  template <class Point>
  index_t localBinsFromPosition(const Point& point) const {
    return grid_helper::getLocalBinIndices(point, m_axes);
  }

  index_t localBinsFromGlobalBin(size_t bin) const {
    return grid_helper::getLocalBinIndices(bin, m_axes);
  }

  template <class Point>
  index_t localBinsFromLowerLeftEdge(const Point& point) const {
    Point shiftedPoint;
    point_t width = grid_helper::getWidth(m_axes);
    for (size_t i = 0; i < DIM; i++) {
      shiftedPoint[i] = point[i] + width[i] / 2;
    }
    return grid_helper::getLocalBinIndices(shiftedPoint, m_axes);
  }

  point_t lowerLeftBinEdge(const index_t& localBins) const {
    return grid_helper::getLowerLeftBinEdge(localBins, m_axes);
  }

  point_t upperRightBinEdge(const index_t& localBins) const {
    return grid_helper::getUpperRightBinEdge(localBins, m_axes);
  }

  index_t numLocalBins() const { return grid_helper::getNBins(m_axes); }

  point_t minPosition() const { return grid_helper::getMin(m_axes); }

  point_t maxPosition() const { return grid_helper::getMax(m_axes); }

  void setExteriorBins(const value_type& value) {
    for (size_t index : grid_helper::exteriorBinIndices(m_axes)) {
      at(index) = value;
    }
  }

  template <
      class Point, typename U = T,
      typename = std::enable_if_t<can_interpolate<
          Point, std::array<double, DIM>, std::array<double, DIM>, U>::value>>
  T interpolate(const Point& point) const {
    // there are 2^DIM corner points used during the interpolation
    constexpr size_t nCorners = 1 << DIM;

    // construct vector of pairs of adjacent bin centers and values
    std::array<value_type, nCorners> neighbors;

    // get local indices for current bin
    // value of bin is interpreted as being the field value at its lower left
    // corner
    const auto& llIndices = localBinsFromPosition(point);

    // get global indices for all surrounding corner points
    const auto& closestIndices = rawClosestPointsIndices(llIndices);

    // get values on grid points
    size_t i = 0;
    for (size_t index : closestIndices) {
      neighbors.at(i++) = at(index);
    }

    return Acts::interpolate(point, lowerLeftBinEdge(llIndices),
                             upperRightBinEdge(llIndices), neighbors);
  }

  template <class Point>
  bool isInside(const Point& position) const {
    return grid_helper::isInside(position, m_axes);
  }

  detail::GlobalNeighborHoodIndices<DIM> neighborHoodIndices(
      const index_t& localBins, size_t size = 1u) const {
    return grid_helper::neighborHoodIndices(localBins, size, m_axes);
  }

  size_t size() const {
    index_t nBinsArray = numLocalBins();
    // add under-and overflow bins for each axis and multiply all bins
    return std::accumulate(
        nBinsArray.begin(), nBinsArray.end(), 1,
        [](const size_t& a, const size_t& b) { return a * (b + 2); });
  }

  std::array<const IAxis*, DIM> axes() const {
    return grid_helper::getAxes(m_axes);
  }

 private:
  std::tuple<Axes...> m_axes;
  std::vector<T> m_values;

  // Part of closestPointsIndices that goes after local bins resolution.
  // Used as an interpolation performance optimization, but not exposed as it
  // doesn't make that much sense from an API design standpoint.
  detail::GlobalNeighborHoodIndices<DIM> rawClosestPointsIndices(
      const index_t& localBins) const {
    return grid_helper::closestPointsIndices(localBins, m_axes);
  }
};
}  // namespace detail

}  // namespace Acts