MDHistoWorkspaceIterator.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidDataObjects/MDHistoWorkspaceIterator.h"
#include "MantidGeometry/MDGeometry/IMDDimension.h"
#include "MantidKernel/System.h"
#include "MantidKernel/Utils.h"
#include "MantidKernel/VMD.h"
#include <algorithm>
#include <boost/math/special_functions/round.hpp>
#include <functional>
#include <numeric>
#include <utility>
 
using namespace Mantid::Kernel;
using namespace Mantid::API;
using Mantid::Geometry::IMDDimension_const_sptr;
 
namespace {
// To get the rounded difference, we need to take into account precision issues
// which arise when
// the bin centres match
Mantid::coord_t getDExact(Mantid::coord_t location, Mantid::coord_t origin, Mantid::coord_t binWidth) {
  auto dExact = (location - origin) / binWidth;
  const auto tolerance = Mantid::coord_t(1e-5);
 
  // Expl. of the steps below
  // -1     -0.5       0      0.5        1   where integer values are bin
  // boundaries and half-values
  //  |       |        |       |         |   are bin centres
  //                               |
  //                              Location
  // Goal: Find distance to closest bin centre:
  // 1. Get the predecimal value of the location: here 0
  // 2. Add and subtract 0.5 to find two possible centres: here -0.5 adn + 0.5
  // 3. Find the centre which minimizes the difference to the location
  // 4. Keep the difference (distance) that was calculation in the step above
 
  // Get the two centre indices which are the closest to dExact.
  const auto centre1 = static_cast<Mantid::coord_t>(size_t(dExact)) + 0.5;
  const auto centre2 = static_cast<Mantid::coord_t>(size_t(dExact)) - 0.5;
 
  // Calculate the distance of the location to the centres
  const auto diff1 = static_cast<Mantid::coord_t>(fabs(centre1 - dExact));
  const auto diff2 = static_cast<Mantid::coord_t>(fabs(centre2 - dExact));
 
  const auto difference = diff1 < diff2 ? diff1 : diff2;
 
  // If the differnce is within the tolerance, ie the location is essentially on
  // the centre, then
  // we want to be on the left of that centre. The correct index for when we are
  // on the centre,
  // is the lower next integer value.
  if (difference < tolerance) {
    if (difference == 0.0f) {
      dExact -= tolerance; // When we have an exact hit in the centre, give it a
                           // nudge to the lower bin boundary
    } else {
      dExact -= 2 * difference; // Need to subtract twice the differnce, in
                                // order to be guaranteed below the centre
    }
  }
  return dExact;
}
} // namespace
 
namespace Mantid::DataObjects {
 
//----------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------
MDHistoWorkspaceIterator::MDHistoWorkspaceIterator(const MDHistoWorkspace_const_sptr &workspace,
                                                   Mantid::Geometry::MDImplicitFunction *function, size_t beginPos,
                                                   size_t endPos)
    : m_skippingPolicy(new SkipMaskedBins(this)) {
  this->init(workspace.get(), function, beginPos, endPos);
}
 
MDHistoWorkspaceIterator::MDHistoWorkspaceIterator(const MDHistoWorkspace *workspace,
                                                   Mantid::Geometry::MDImplicitFunction *function, size_t beginPos,
                                                   size_t endPos)
    : m_skippingPolicy(new SkipMaskedBins(this)) {
  this->init(workspace, function, beginPos, endPos);
}
 
MDHistoWorkspaceIterator::MDHistoWorkspaceIterator(const MDHistoWorkspace_const_sptr &workspace,
                                                   SkippingPolicy *skippingPolicy,
                                                   Mantid::Geometry::MDImplicitFunction *function, size_t beginPos,
                                                   size_t endPos)
    : m_skippingPolicy(skippingPolicy) {
  this->init(workspace.get(), function, beginPos, endPos);
}
 
//----------------------------------------------------------------------------------------------
MDHistoWorkspaceIterator::MDHistoWorkspaceIterator(const MDHistoWorkspace *workspace, SkippingPolicy *skippingPolicy,
                                                   Mantid::Geometry::MDImplicitFunction *function, size_t beginPos,
                                                   size_t endPos)
    : m_skippingPolicy(skippingPolicy) {
  this->init(workspace, function, beginPos, endPos);
}
 
void MDHistoWorkspaceIterator::init(const MDHistoWorkspace *workspace, Mantid::Geometry::MDImplicitFunction *function,
                                    size_t beginPos, size_t endPos) {
  m_ws = workspace;
  if (m_ws == nullptr)
    throw std::invalid_argument("MDHistoWorkspaceIterator::ctor(): NULL workspace given.");
 
  m_begin = beginPos;
  m_pos = m_begin;
  m_function.reset(function);
 
  m_max = endPos;
  if (m_max > m_ws->getNPoints())
    m_max = m_ws->getNPoints();
  if (m_max < m_pos)
    throw std::invalid_argument("MDHistoWorkspaceIterator::ctor(): End point "
                                "given is before the start point.");
 
  m_nd = m_ws->getNumDims();
  m_center = new coord_t[m_nd];
  m_origin = new coord_t[m_nd];
  m_binWidth = new coord_t[m_nd];
  m_index = new size_t[m_nd];
  m_indexMax = new size_t[m_nd];
  m_indexMaker = new size_t[m_nd];
  Utils::NestedForLoop::SetUp(m_nd, m_index, 0);
  // Initalize all these values
  for (size_t d = 0; d < m_nd; d++) {
    IMDDimension_const_sptr dim = m_ws->getDimension(d);
    m_center[d] = 0;
    m_origin[d] = dim->getMinimum();
    m_binWidth[d] = dim->getBinWidth();
    m_indexMax[d] = dim->getNBins();
  }
  Utils::NestedForLoop::SetUpIndexMaker(m_nd, m_indexMaker, m_indexMax);
 
  // Initialize the current index from the start position.
  Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker, m_indexMax, m_index);
 
  // Make sure that the first iteration is at a point inside the implicit
  // function
  if (m_function) {
    // Calculate the center of the 0-th bin
    for (size_t d = 0; d < m_nd; d++)
      m_center[d] = m_origin[d] + 0.5f * m_binWidth[d];
    // Skip on if the first point is NOT contained
    if (!m_function->isPointContained(m_center)) {
      bool didNext = next();
      if (!didNext && this->valid()) {
        throw std::runtime_error("Inconsistency found initializing "
                                 "MDHistoWorkspace iterator: this iterator should be valid, but "
                                 "when tried to skip the "
                                 "first point (not contained) could not iterate to "
                                 "next point.");
      }
    }
  }
 
  // --- Calculate index permutations for neighbour finding face touching ---
  auto temp = std::vector<int64_t>(2 * m_nd);
  m_permutationsFaceTouching.swap(temp);
 
  int64_t offset = 1;
 
  m_permutationsFaceTouching[0] = -1;
  m_permutationsFaceTouching[1] = 1;
 
  // Figure out what possible indexes deltas to generate indexes that are next
  // to the current one.
  for (size_t j = 1; j < m_nd; ++j) {
    offset = offset * static_cast<int64_t>(m_ws->getDimension(j - 1)->getNBins());
 
    m_permutationsFaceTouching[j * 2] = offset;
    m_permutationsFaceTouching[(j * 2) + 1] = -offset;
  }
}
 
//----------------------------------------------------------------------------------------------
MDHistoWorkspaceIterator::~MDHistoWorkspaceIterator() {
  delete[] m_center;
  delete[] m_origin;
  delete[] m_binWidth;
  delete[] m_index;
  delete[] m_indexMax;
  delete[] m_indexMaker;
}
//----------------------------------------------------------------------------------------------
size_t MDHistoWorkspaceIterator::getDataSize() const { return size_t(m_max - m_begin); }
 
//----------------------------------------------------------------------------------------------
void MDHistoWorkspaceIterator::jumpTo(size_t index) { m_pos = uint64_t(index + m_begin); }
 
Mantid::coord_t MDHistoWorkspaceIterator::jumpToNearest(const VMD &fromLocation) {
  std::vector<size_t> indexes(m_nd);
  coord_t sqDiff = 0;
  for (size_t d = 0; d < m_nd; ++d) {
    coord_t dExact = getDExact(fromLocation[d], m_origin[d], m_binWidth[d]);
    size_t dRound = std::lround(dExact); // Round to nearest bin edge.
    if (dRound >= m_indexMax[d]) {
      dRound = m_indexMax[d] - 1;
    }
    sqDiff += (dExact - coord_t(dRound)) * (dExact - coord_t(dRound)) * m_binWidth[d] * m_binWidth[d];
    indexes[d] = dRound;
  }
 
  const size_t linearIndex = Utils::NestedForLoop::GetLinearIndex(m_nd, &indexes[0], m_indexMaker);
  this->jumpTo(linearIndex);
  return std::sqrt(sqDiff);
}
 
//----------------------------------------------------------------------------------------------
bool MDHistoWorkspaceIterator::valid() const { return (m_pos < m_max); }
 
//----------------------------------------------------------------------------------------------
bool MDHistoWorkspaceIterator::next() {
  if (m_function) {
    bool allIncremented;
 
    do {
      m_pos++;
      allIncremented = Utils::NestedForLoop::Increment(m_nd, m_index, m_indexMax);
      // Calculate the center
      for (size_t d = 0; d < m_nd; d++) {
        m_center[d] = m_origin[d] + (coord_t(m_index[d]) + 0.5f) * m_binWidth[d];
      }
      // Keep incrementing until you are in the implicit function and not masked
    } while (m_pos < m_max &&
             ((!allIncremented && !m_function->isPointContained(m_center)) || m_skippingPolicy->keepGoing()));
  } else {
    // Keep moving to next position if the current position is masked and
    // still valid.
    do {
      m_pos++;
    } while (m_pos < m_max && m_skippingPolicy->keepGoing());
  }
 
  bool ret = m_pos < m_max;
  // Go through every point;
  return ret;
}
 
//----------------------------------------------------------------------------------------------
bool MDHistoWorkspaceIterator::next(size_t skip) {
  m_pos += skip;
 
  return (m_pos < m_max);
}
 
//----------------------------------------------------------------------------------------------
signal_t MDHistoWorkspaceIterator::getNormalizedSignal() const {
  // What is our normalization factor?
  switch (m_normalization) {
  case NoNormalization:
    return m_ws->getSignalAt(m_pos);
  case VolumeNormalization:
    return m_ws->getSignalAt(m_pos) * m_ws->getInverseVolume();
  case NumEventsNormalization:
    return m_ws->getSignalAt(m_pos) / m_ws->getNumEventsAt(m_pos);
  }
  // Should not reach here
  return std::numeric_limits<signal_t>::quiet_NaN();
}
 
//----------------------------------------------------------------------------------------------
signal_t MDHistoWorkspaceIterator::getNormalizedError() const {
  // What is our normalization factor?
  switch (m_normalization) {
  case NoNormalization:
    return m_ws->getErrorAt(m_pos);
  case VolumeNormalization:
    return m_ws->getErrorAt(m_pos) * m_ws->getInverseVolume();
  case NumEventsNormalization:
    return m_ws->getErrorAt(m_pos) / m_ws->getNumEventsAt(m_pos);
  }
  // Should not reach here
  return std::numeric_limits<signal_t>::quiet_NaN();
}
 
//----------------------------------------------------------------------------------------------
signal_t MDHistoWorkspaceIterator::getSignal() const { return m_ws->getSignalAt(m_pos); }
 
signal_t MDHistoWorkspaceIterator::getError() const { return m_ws->getErrorAt(m_pos); }
//----------------------------------------------------------------------------------------------
std::unique_ptr<coord_t[]> MDHistoWorkspaceIterator::getVertexesArray(size_t &numVertices) const {
  // The MDHistoWorkspace takes care of this
  return m_ws->getVertexesArray(m_pos, numVertices);
}
 
std::unique_ptr<coord_t[]> MDHistoWorkspaceIterator::getVertexesArray(size_t &numVertices, const size_t outDimensions,
                                                                      const bool *maskDim) const {
  // Do the same thing as is done in the MDBoxBase
  UNUSED_ARG(numVertices);
  UNUSED_ARG(outDimensions);
  UNUSED_ARG(maskDim);
  throw std::runtime_error("Not Implemented At present time");
}
 
//----------------------------------------------------------------------------------------------
Mantid::Kernel::VMD MDHistoWorkspaceIterator::getCenter() const {
  // Get the indices
  Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker, m_indexMax, m_index);
  // Find the center
  for (size_t d = 0; d < m_nd; ++d) {
    m_center[d] = m_origin[d] + (coord_t(m_index[d]) + 0.5f) * m_binWidth[d];
  }
  return VMD(m_nd, m_center);
}
 
VecMDExtents MDHistoWorkspaceIterator::getBoxExtents() const {
 
  // Get the indexes.
  Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker, m_indexMax, m_index);
  VecMDExtents extents(m_nd);
  // Find the extents.
  for (size_t d = 0; d < m_nd; ++d) {
    const coord_t min = m_origin[d] + coord_t(m_index[d]) * m_binWidth[d]; // Min in d
    const coord_t max = min + m_binWidth[d];                               // Max in d
    extents[d] = MDExtentPair(min, max);
  }
 
  return extents;
}
 
//----------------------------------------------------------------------------------------------
size_t MDHistoWorkspaceIterator::getNumEvents() const { return static_cast<size_t>(this->getNumEventsFraction()); }
 
//----------------------------------------------------------------------------------------------
signal_t MDHistoWorkspaceIterator::getNumEventsFraction() const { return m_ws->getNumEventsAt(m_pos); }
 
//----------------------------------------------------------------------------------------------
uint16_t MDHistoWorkspaceIterator::getInnerExpInfoIndex(size_t /*index*/) const {
  return 0;
  // throw std::runtime_error("MDHistoWorkspaceIterator: No events are
  // contained, so it is not possible to return inner associated experiment-info index.");
}
 
uint16_t MDHistoWorkspaceIterator::getInnerGoniometerIndex(size_t /*index*/) const {
  return 0;
  // throw std::runtime_error("MDHistoWorkspaceIterator: No events are
  // contained, so it is not possible to return inner goniometer index.");
}
 
int32_t MDHistoWorkspaceIterator::getInnerDetectorID(size_t /*index*/) const {
  return 0;
  // throw std::runtime_error("MDHistoWorkspaceIterator: No events are
  // contained, so it is not possible to return inner detector ID.");
}
 
coord_t MDHistoWorkspaceIterator::getInnerPosition(size_t /*index*/, size_t dimension) const {
  return m_ws->getCenter(m_pos)[dimension];
}
 
signal_t MDHistoWorkspaceIterator::getInnerSignal(size_t /*index*/) const { return m_ws->getSignalAt(m_pos); }
 
signal_t MDHistoWorkspaceIterator::getInnerError(size_t /*index*/) const { return m_ws->getErrorAt(m_pos); }
 
bool MDHistoWorkspaceIterator::getIsMasked() const { return m_ws->getIsMaskedAt(m_pos); }
 
size_t MDHistoWorkspaceIterator::getLinearIndex() const { return m_pos; }
 
std::vector<size_t> MDHistoWorkspaceIterator::findNeighbourIndexes() const {
 
  return this->findNeighbourIndexesByWidth(3 /*immediate neighbours only*/);
}
 
std::vector<size_t> MDHistoWorkspaceIterator::findNeighbourIndexesFaceTouching() const {
  Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker, m_indexMax, m_index);
 
  std::vector<size_t> neighbourIndexes; // Accumulate neighbour indexes.
  std::vector<int> widths(m_nd, 3);     // Face touching width is always 3 in each dimension
  for (auto permutation : m_permutationsFaceTouching) {
    if (permutation == 0) {
      continue;
    }
 
    size_t neighbour_index = m_pos + permutation;
    if (neighbour_index < m_ws->getNPoints() &&
        Utils::isNeighbourOfSubject(m_nd, neighbour_index, m_index, m_indexMaker, m_indexMax, widths)) {
      neighbourIndexes.emplace_back(neighbour_index);
    }
  }
  return neighbourIndexes;
}
 
bool MDHistoWorkspaceIterator::isWithinBounds(size_t index) const { return index >= m_begin && index < m_max; }
 
std::vector<int64_t> MDHistoWorkspaceIterator::createPermutations(const std::vector<int> &widths) const {
  // look-up
  auto it = m_permutationsVertexTouchingMap.find(widths);
  if (it == m_permutationsVertexTouchingMap.end()) {
 
    if (widths[0] % 2 == 0) {
      throw std::invalid_argument("MDHistoWorkspaceIterator::"
                                  "findNeighbourIndexesByWidth, width must "
                                  "always be an odd number");
    }
    if (widths.size() != m_nd) {
      throw std::invalid_argument("MDHistoWorkspaceIterator::"
                                  "findNeighbourIndexesByWidth, size of widths "
                                  "must be the same as the number of "
                                  "dimensions.");
    }
 
    int64_t offset = 1;
 
    // Size of block will be width ^ nd
    std::vector<int64_t> permutationsVertexTouching;
    // Calculate maximum permutations size.
    int product = std::accumulate(widths.begin(), widths.end(), 1, std::multiplies<int>());
    permutationsVertexTouching.reserve(product);
 
    int centreIndex = widths[0] / 2; // Deliberately truncate to get centre index
 
    for (int i = 0; i < widths[0]; ++i) {
      // for width = 3 : -1, 0, 1
      // for width = 5 : -2, -1, 0, 1, 2
      permutationsVertexTouching.emplace_back(centreIndex - i);
    }
 
    // Figure out what possible indexes deltas to generate indexes that are next
    // to the current one.
    for (size_t j = 1; j < m_nd; ++j) {
      offset = offset * static_cast<int64_t>(m_ws->getDimension(j - 1)->getNBins());
 
      size_t nEntries = permutationsVertexTouching.size();
      for (int k = 1; k <= widths[j] / 2; ++k) {
        for (size_t m = 0; m < nEntries; m++) {
          permutationsVertexTouching.emplace_back((offset * k) + permutationsVertexTouching[m]);
          permutationsVertexTouching.emplace_back((offset * k * (-1)) + permutationsVertexTouching[m]);
        }
      }
    }
 
    m_permutationsVertexTouchingMap.insert(std::make_pair(widths, permutationsVertexTouching));
  }
 
  // In either case, get the result.
  return m_permutationsVertexTouchingMap[widths];
}
 
std::vector<size_t> MDHistoWorkspaceIterator::findNeighbourIndexesByWidth(const int &width) const {
 
  return this->findNeighbourIndexesByWidth(std::vector<int>(m_nd, width));
}
 
std::vector<size_t> MDHistoWorkspaceIterator::findNeighbourIndexesByWidth(const std::vector<int> &widths) const {
 
  // Find existing or create required index permutations.
  std::vector<int64_t> permutationsVertexTouching = createPermutations(widths);
 
  Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker, m_indexMax, m_index);
 
  // Filter out indexes that are are not actually neighbours.
  // Accumulate neighbour indexes.
  std::vector<size_t> neighbourIndexes(permutationsVertexTouching.size());
  size_t nextFree = 0;
  for (auto permutation : permutationsVertexTouching) {
    if (permutation == 0) {
      continue;
    }
 
    size_t neighbour_index = m_pos + permutation;
    if (neighbour_index < m_ws->getNPoints() &&
        Utils::isNeighbourOfSubject(m_nd, neighbour_index, m_index, m_indexMaker, m_indexMax, widths)) {
      neighbourIndexes[nextFree++] = neighbour_index;
    }
  }
  neighbourIndexes.resize(nextFree);
 
  // Remove duplicates
  std::sort(neighbourIndexes.begin(), neighbourIndexes.end());
  neighbourIndexes.erase(std::unique(neighbourIndexes.begin(), neighbourIndexes.end()), neighbourIndexes.end());
  return neighbourIndexes;
}
 
std::pair<std::vector<size_t>, std::vector<bool>>
MDHistoWorkspaceIterator::findNeighbourIndexesByWidth1D(const int &width, const int &width_dimension) const {
 
  std::vector<int> widths;
  for (size_t dimension = 0; dimension < m_nd; ++dimension) {
    if (static_cast<int>(dimension) == width_dimension) {
      widths.emplace_back(width);
    } else {
      widths.emplace_back(1);
    }
  }
 
  // Find existing or create required index permutations.
  std::vector<int64_t> permutationsVertexTouching = createPermutations(widths);
 
  std::vector<bool> indexValidity(permutationsVertexTouching.size(), false);
 
  Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker, m_indexMax, m_index);
 
  std::sort(permutationsVertexTouching.begin(), permutationsVertexTouching.end());
 
  // Accumulate neighbour indexes.
  // Record indexes as valid only if they are actually neighbours.
  std::vector<size_t> neighbourIndexes(permutationsVertexTouching.size());
  for (size_t i = 0; i < permutationsVertexTouching.size(); ++i) {
 
    size_t neighbour_index = m_pos + permutationsVertexTouching[i];
    neighbourIndexes[i] = neighbour_index;
    if (neighbour_index < m_ws->getNPoints() &&
        Utils::isNeighbourOfSubject(m_nd, neighbour_index, m_index, m_indexMaker, m_indexMax, widths)) {
      indexValidity[i] = true;
    }
  }
 
  return std::make_pair(neighbourIndexes, indexValidity);
}
 
size_t MDHistoWorkspaceIterator::permutationCacheSize() const { return m_permutationsVertexTouchingMap.size(); }
 
} // namespace Mantid::DataObjects