MDGridBox.hxx

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/CoordTransformDistance.h"
#include "MantidDataObjects/MDBox.h"
#include "MantidDataObjects/MDEvent.h"
#include "MantidDataObjects/MDGridBox.h"
#include "MantidKernel/FunctionTask.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/Task.h"
#include "MantidKernel/ThreadPool.h"
#include "MantidKernel/ThreadScheduler.h"
#include "MantidKernel/ThreadSchedulerMutexes.h"
#include "MantidKernel/Timer.h"
#include "MantidKernel/Utils.h"
#include "MantidKernel/WarningSuppressions.h"
#include <boost/math/special_functions/round.hpp>
#include <optional>
#include <ostream>
 
// These pragmas ignores the warning in the ctor where "d<nd-1" for nd=1.
// This is okay (though would be better if it were for only that function
#if (defined(__INTEL_COMPILER))
#pragma warning disable 186
#elif defined(__GNUC__)
#pragma GCC diagnostic ignored "-Wtype-limits"
#endif
 
namespace Mantid {
namespace DataObjects {
 
 
//-----------------------------------------------------------------------------------------------
TMDE(MDGridBox)::MDGridBox(API::BoxController *const bc, const uint32_t depth,
                           const std::vector<Mantid::Geometry::MDDimensionExtents<coord_t>> &extentsVector)
    : MDBoxBase<MDE, nd>(bc, depth, UNDEF_SIZET, extentsVector), numBoxes(0), m_Children(), diagonalSquared(0.f),
      nPoints(0) {
  initGridBox();
}
 
TMDE(MDGridBox)::MDGridBox(std::shared_ptr<API::BoxController> &bc, const uint32_t depth,
                           const std::vector<Mantid::Geometry::MDDimensionExtents<coord_t>> &extentsVector)
    : MDBoxBase<MDE, nd>(bc.get(), depth, UNDEF_SIZET, extentsVector), numBoxes(0), m_Children(), diagonalSquared(0.f),
      nPoints(0) {
  initGridBox();
}
 
template <typename MDE, size_t nd> size_t MDGridBox<MDE, nd>::initGridBox() {
  if (!this->m_BoxController)
    throw std::runtime_error("MDGridBox::ctor(): No BoxController specified in box.");
 
  // How many is it split?
  // If we are at the top level and we have a specific top level split, then set
  // it.
  std::optional<std::vector<size_t>> splitTopInto = this->m_BoxController->getSplitTopInto();
  if (this->getDepth() == 0 && splitTopInto) {
    for (size_t d = 0; d < nd; d++)
      split[d] = splitTopInto.value()[d];
  } else {
    for (size_t d = 0; d < nd; d++)
      split[d] = this->m_BoxController->getSplitInto(d);
  }
 
  // Compute sizes etc.
  size_t tot = computeSizesFromSplit();
  if (tot == 0)
    throw std::runtime_error("MDGridBox::ctor(): Invalid splitting criterion (one was zero).");
  return tot;
}
 
//-----------------------------------------------------------------------------------------------
TMDE(MDGridBox)::MDGridBox(MDBox<MDE, nd> *box)
    : MDBoxBase<MDE, nd>(*box, box->getBoxController()), split(), splitCumul(), m_SubBoxSize(), numBoxes(0),
      m_Children(), diagonalSquared(0.f), nPoints(0) {
  size_t totalSize = initGridBox();
 
  double ChildVol(1);
  for (size_t d = 0; d < nd; d++)
    ChildVol *= m_SubBoxSize[d];
 
  // Splitting an input MDBox requires creating a bunch of children
  fillBoxShell(totalSize, coord_t(1. / ChildVol));
 
  // Prepare to distribute the events that were in the box before, this will
  // load missing events from HDD in file based ws if there are some.
  const std::vector<MDE> &events = box->getConstEvents();
  // just add event to the existing internal box
  for (const auto &evnt : events)
    addEvent(evnt);
 
  // Copy the cached numbers from the incoming box. This is quick - don't need
  // to refresh cache
  this->nPoints = box->getNPoints();
 
  // Clear the old box and delete it from disk buffer if one is used.
  box->clear();
}
template <typename MDE, size_t nd>
void MDGridBox<MDE, nd>::fillBoxShell(const size_t tot, const coord_t ChildInverseVolume) {
  // Create the array of MDBox contents.
  this->m_Children.clear();
  this->m_Children.reserve(tot);
  this->numBoxes = tot;
 
  size_t indices[nd];
  for (size_t d = 0; d < nd; d++)
    indices[d] = 0;
 
  // get inital free ID for the boxes, which would be created by this command
  // Splitting an input MDBox requires creating a bunch of children
  // But the IDs of these children MUST be sequential. Hence the critical block
  // within claimIDRange,
  // which would produce sequental ranges in multithreaded environment
  size_t ID0 = this->m_BoxController->claimIDRange(tot);
 
  for (size_t i = 0; i < tot; i++) {
    // Create the box
    // (Increase the depth of this box to one more than the parent (this))
    auto splitBox = new MDBox<MDE, nd>(this->m_BoxController, this->m_depth + 1, UNDEF_SIZET, size_t(ID0 + i));
    // This MDGridBox is the parent of the new child.
    splitBox->setParent(this);
 
    // Set the extents of this box.
    for (size_t d = 0; d < nd; d++) {
      double min = double(this->extents[d].getMin()) + double(indices[d]) * m_SubBoxSize[d];
      double max = min + m_SubBoxSize[d];
      splitBox->setExtents(d, min, max);
    }
    splitBox->setInverseVolume(ChildInverseVolume); // Set the cached inverse volume
    m_Children.emplace_back(splitBox);
 
    // Increment the indices, rolling back as needed
    indices[0]++;
    for (size_t d = 0; d < nd - 1; d++) // This is not run if nd=1; that's okay,
                                        // you can ignore the warning
    {
      if (indices[d] >= split[d]) {
        indices[d] = 0;
        indices[d + 1]++;
      }
    }
  } // for each box
}
 
//-----------------------------------------------------------------------------------------------
TMDE(MDGridBox)::MDGridBox(const MDGridBox<MDE, nd> &other, Mantid::API::BoxController *const otherBC)
    : MDBoxBase<MDE, nd>(other, otherBC), numBoxes(other.numBoxes), m_Children(),
      diagonalSquared(other.diagonalSquared), nPoints(other.nPoints) {
  for (size_t d = 0; d < nd; d++) {
    split[d] = other.split[d];
    splitCumul[d] = other.splitCumul[d];
    m_SubBoxSize[d] = other.m_SubBoxSize[d];
    m_SubBoxSizeInv[d] = other.m_SubBoxSizeInv[d];
  }
  // Copy all the boxes
  m_Children.clear();
  m_Children.reserve(numBoxes);
  for (size_t i = 0; i < other.m_Children.size(); i++) {
    API::IMDNode *otherBox = other.m_Children[i];
    const MDBox<MDE, nd> *otherMDBox = dynamic_cast<const MDBox<MDE, nd> *>(otherBox);
    const MDGridBox<MDE, nd> *otherMDGridBox = dynamic_cast<const MDGridBox<MDE, nd> *>(otherBox);
    if (otherMDBox) {
      auto newBox = new MDBox<MDE, nd>(*otherMDBox, otherBC);
      newBox->setParent(this);
      m_Children.emplace_back(newBox);
    } else if (otherMDGridBox) {
      auto newBox = new MDGridBox<MDE, nd>(*otherMDGridBox, otherBC);
      newBox->setParent(this);
      m_Children.emplace_back(newBox);
    } else {
      throw std::runtime_error("MDGridBox::copy_ctor(): an unexpected child box type was found.");
    }
  }
}
 
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::transformDimensions(std::vector<double> &scaling, std::vector<double> &offset) {
  MDBoxBase<MDE, nd>::transformDimensions(scaling, offset);
  this->computeSizesFromSplit();
}
 
//-----------------------------------------------------------------------------------------------
TMDE(size_t MDGridBox)::computeSizesFromSplit() {
  // Do some computation based on how many splits per each dim.
  size_t tot = 1;
  double diagSum(0);
  for (size_t d = 0; d < nd; d++) {
    // Cumulative multiplier, for indexing
    splitCumul[d] = tot;
    tot *= split[d];
    // Length of the side of a box in this dimension
    m_SubBoxSize[d] = static_cast<double>(this->extents[d].getSize()) / static_cast<double>(split[d]);
    m_SubBoxSizeInv[d] = 1.0 / m_SubBoxSize[d];
    // Accumulate the squared diagonal length.
    diagSum += m_SubBoxSize[d] * m_SubBoxSize[d];
  }
  diagonalSquared = static_cast<coord_t>(diagSum);
 
  return tot;
}
 
//-----------------------------------------------------------------------------------------------
TMDE(MDGridBox)::~MDGridBox() {
  // Delete all contained boxes (this should fire the MDGridBox destructors
  // recursively).
  auto it = m_Children.begin();
  for (; it != m_Children.end(); ++it)
    delete *it;
  m_Children.clear();
}
 
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::clear() {
  this->m_signal = 0.0;
  this->m_errorSquared = 0.0;
  auto it = m_Children.begin();
  for (; it != m_Children.end(); ++it) {
    (*it)->clear();
  }
}
 
//-----------------------------------------------------------------------------------------------
TMDE(size_t MDGridBox)::getNumDims() const { return nd; }
 
//-----------------------------------------------------------------------------------------------
TMDE(size_t MDGridBox)::getDataInMemorySize() const {
  size_t nPoints(0);
  for (size_t i = 0; i < numBoxes; i++)
    nPoints += m_Children[i]->getDataInMemorySize();
  return nPoints;
}
 
//-----------------------------------------------------------------------------------------------
TMDE(size_t MDGridBox)::getNumMDBoxes() const {
  size_t total = 0;
  auto it = m_Children.begin();
  for (; it != m_Children.end(); ++it) {
    total += (*it)->getNumMDBoxes();
  }
  return total;
}
 
//-----------------------------------------------------------------------------------------------
TMDE(size_t MDGridBox)::getNumChildren() const { return numBoxes; }
 
//-----------------------------------------------------------------------------------------------
TMDE(API::IMDNode *MDGridBox)::getChild(size_t index) { return m_Children[index]; }
 
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::setChildren(const std::vector<API::IMDNode *> &otherBoxes, const size_t indexStart,
                                  const size_t indexEnd) {
  m_Children.clear();
  m_Children.reserve(indexEnd - indexStart + 1);
  auto it = otherBoxes.begin() + indexStart;
  auto it_end = otherBoxes.begin() + indexEnd;
  // Set the parent of each new child box.
  for (; it != it_end; it++) {
    m_Children.emplace_back(dynamic_cast<MDBoxBase<MDE, nd> *>(*it));
    m_Children.back()->setParent(this);
  }
  numBoxes = m_Children.size();
}
 
//-----------------------------------------------------------------------------------------------
TMDE(inline size_t MDGridBox)::getLinearIndex(size_t *indices) const {
  size_t out_linear_index = 0;
  for (size_t d = 0; d < nd; d++)
    out_linear_index += (indices[d] * splitCumul[d]);
  return out_linear_index;
}
 
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::refreshCache(Kernel::ThreadScheduler *ts) {
  // Clear your total
  nPoints = 0;
  this->m_signal = 0;
  this->m_errorSquared = 0;
  this->m_totalWeight = 0;
 
  if (!ts) {
    //--------- Serial -----------
    for (MDBoxBase<MDE, nd> *ibox : m_Children) {
 
      // Refresh the cache (does nothing for MDBox)
      ibox->refreshCache();
 
      // Add up what's in there
      nPoints += ibox->getNPoints();
      this->m_signal += ibox->getSignal();
      this->m_errorSquared += ibox->getErrorSquared();
      this->m_totalWeight += ibox->getTotalWeight();
    }
  } else {
    //---------- Parallel refresh --------------
    throw std::runtime_error("Not implemented");
  }
}
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::calculateGridCaches() {
  // Clear your total
  nPoints = 0;
  this->m_signal = 0;
  this->m_errorSquared = 0;
  this->m_totalWeight = 0;
 
  for (MDBoxBase<MDE, nd> *ibox : m_Children) {
 
    // does nothing for MDBox
    ibox->calculateGridCaches();
 
    // Add up what's in there
    nPoints += ibox->getNPoints();
    this->m_signal += ibox->getSignal();
    this->m_errorSquared += ibox->getErrorSquared();
    this->m_totalWeight += ibox->getTotalWeight();
  }
}
//-----------------------------------------------------------------------------------------------
TMDE(std::vector<MDE> *MDGridBox)::getEventsCopy() {
  auto out = new std::vector<MDE>();
  // Make the copy
  // out->insert(out->begin(), data.begin(), data.end());
  return out;
}
 
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::getBoxes(std::vector<API::IMDNode *> &outBoxes, size_t maxDepth, bool leafOnly) {
  // Add this box, unless we only want the leaves
  if (!leafOnly)
    outBoxes.emplace_back(this);
 
  if (this->getDepth() < maxDepth) {
    for (API::IMDNode *child : m_Children) {
      // Recursively go deeper, if needed
      child->getBoxes(outBoxes, maxDepth, leafOnly);
    }
  } else {
    // Oh, we reached the max depth and want only leaves.
    // ... so we consider this box to be a leaf too.
    if (leafOnly)
      outBoxes.emplace_back(this);
  }
}
 
GNU_DIAG_OFF("maybe-uninitialized")
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::getBoxes(std::vector<API::IMDNode *> &outBoxes, size_t maxDepth, bool leafOnly,
                               Mantid::Geometry::MDImplicitFunction *function) {
  // Add this box, unless we only want the leaves
  if (!leafOnly)
    outBoxes.emplace_back(this);
 
  if (this->getDepth() < maxDepth) {
    // OK, let's look for children that are either touching or completely
    // contained by the implicit function.
 
    // The number of vertices in each dimension is the # split[d] + 1
    size_t vertices_max[nd];
    Kernel::Utils::NestedForLoop::SetUp(nd, vertices_max, 0);
 
    // Total number of vertices for all the boxes
    size_t numVertices = 1;
    for (size_t d = 0; d < nd; ++d) {
      vertices_max[d] = split[d] + 1;
      numVertices *= vertices_max[d];
    }
 
    // The function is limited by this many planes
    size_t numPlanes = function->getNumPlanes();
 
    // This array will hold whether each vertex is contained by each plane.
    auto vertexContained = new bool[numVertices * numPlanes];
 
    // The index to the vertex in each dimension
    size_t vertexIndex[nd];
    Kernel::Utils::NestedForLoop::SetUp(nd, vertexIndex, 0);
    // To get indexes in the array of vertexes
    size_t vertexIndexMaker[nd];
    Kernel::Utils::NestedForLoop::SetUpIndexMaker(nd, vertexIndexMaker, vertices_max);
    // To get indexes in the array of BOXES
    size_t boxIndexMaker[nd];
    Kernel::Utils::NestedForLoop::SetUpIndexMaker(nd, boxIndexMaker, split);
 
    size_t linearVertexIndex = 0;
    for (linearVertexIndex = 0; linearVertexIndex < numVertices; linearVertexIndex++) {
      // Get the nd-dimensional index
      Kernel::Utils::NestedForLoop::GetIndicesFromLinearIndex(nd, linearVertexIndex, vertexIndexMaker, vertices_max,
                                                              vertexIndex);
 
      // Coordinates of this vertex
      coord_t vertexCoord[nd];
      for (size_t d = 0; d < nd; ++d)
        vertexCoord[d] = this->extents[d].getMin() + coord_t(double(vertexIndex[d]) * m_SubBoxSize[d]);
 
      // Now check each plane to see if the vertex is bounded by it
      for (size_t p = 0; p < numPlanes; p++) {
        // Save whether this vertex is contained by this plane
        vertexContained[p * numVertices + linearVertexIndex] = function->getPlane(p).isPointInside(vertexCoord);
      }
    }
 
    // OK, now we have an array saying which vertex is contained by which plane.
 
    // This is the number of vertices for each box, e.g. 8 in 3D
    size_t verticesPerBox = 1 << nd;
 
    /* There is a fixed relationship betwen a vertex (in a linear index) and its
     * neighbors for a given box. This array calculates this:  */
    auto vertexNeighborsOffsets = new size_t[verticesPerBox];
 
    for (size_t i = 0; i < verticesPerBox; i++) {
      // Index (in n-dimensions) of this neighbor)
      size_t vertIndex[nd];
      for (size_t d = 0; d < nd; d++) {
        vertIndex[d] = 0;
        // Use a bit mask to iterate through the 2^nd neighbor options
        size_t mask = size_t{1} << d;
        if (i & mask)
          vertIndex[d] = 1;
      }
      size_t linIndex = Kernel::Utils::NestedForLoop::GetLinearIndex(nd, vertIndex, vertexIndexMaker);
      vertexNeighborsOffsets[i] = linIndex;
    }
 
    // Go through all the boxes
    size_t boxIndex[nd];
    Kernel::Utils::NestedForLoop::SetUp(nd, boxIndex, 0);
 
    bool allDone = false;
    while (!allDone) {
      // Find the linear index into the BOXES array.
      size_t boxLinearIndex = Kernel::Utils::NestedForLoop::GetLinearIndex(nd, boxIndex, boxIndexMaker);
      API::IMDNode *box = m_Children[boxLinearIndex];
 
      //        std::cout << "Box at " << Strings::join(boxIndex, boxIndex+nd,
      //        ", ")
      //              << " (" << box->getExtentsStr() << ") ";
 
      // Find the linear index of the upper left vertex of the box.
      // (note that we're using the VERTEX index maker to find the linear index
      // in that LARGER array)
      size_t vertLinearIndex = Kernel::Utils::NestedForLoop::GetLinearIndex(nd, boxIndex, vertexIndexMaker);
 
      // OK, now its time to see if the box is touching or contained or out of
      // it.
      // Recall that:
      //  - if a plane has NO vertices, then the box DOES NOT TOUCH
      //  - if EVERY plane has EVERY vertex, then the box is CONTAINED
      //  - if EVERY plane has at least one vertex, then the box is TOUCHING
 
      size_t numPlanesWithAllVertexes = 0;
 
      bool boxIsNotTouching = false;
 
      // Go plane by plane
      for (size_t p = 0; p < numPlanes; p++) {
        size_t numVertexesInThisPlane = 0;
        // Evaluate the 2^nd vertexes for this box.
        for (size_t i = 0; i < verticesPerBox; i++) {
          // (the index of the vertex is) = vertLinearIndex +
          // vertexNeighborsOffsets[i]
          if (vertexContained[p * numVertices + vertLinearIndex + vertexNeighborsOffsets[i]])
            numVertexesInThisPlane++;
        }
 
        // Plane with no vertexes = NOT TOUCHING. You can exit now
        if (numVertexesInThisPlane == 0) {
          boxIsNotTouching = true;
          break;
        }
 
        // Plane has all the vertexes
        if (numVertexesInThisPlane == verticesPerBox)
          numPlanesWithAllVertexes++;
      } // (for each plane)
 
      // Is there a chance that the box is contained?
      if (!boxIsNotTouching) {
 
        if (numPlanesWithAllVertexes == numPlanes) {
          // All planes have all vertexes
          // The box is FULLY CONTAINED
          // So we can get ALL children and don't need to check the implicit
          // function
          box->getBoxes(outBoxes, maxDepth, leafOnly);
        } else {
          // There is a chance the box is touching. Keep checking with implicit
          // functions
          box->getBoxes(outBoxes, maxDepth, leafOnly, function);
        }
      } else {
        //          std::cout << " is not touching at all.\n";
      }
 
      // Move on to the next box in the list
      allDone = Kernel::Utils::NestedForLoop::Increment(nd, boxIndex, split);
    }
 
    // Clean up.
    delete[] vertexContained;
    delete[] vertexNeighborsOffsets;
 
  } // Not at max depth
  else {
    // Oh, we reached the max depth and want only leaves.
    // ... so we consider this box to be a leaf too.
    if (leafOnly)
      outBoxes.emplace_back(this);
  }
}
GNU_DIAG_ON("maybe-uninitialized")
//-----------------------------------------------------------------------------------------------
template <typename MDE, size_t nd> const API::IMDNode *MDGridBox<MDE, nd>::getBoxAtCoord(const coord_t *coords) {
  size_t index = 0;
  for (size_t d = 0; d < nd; d++) {
    coord_t x = coords[d];
    int i = int((x - this->extents[d].getMin()) / m_SubBoxSize[d]);
    // NOTE: No bounds checking is done (for performance).
    // Accumulate the index
    index += (i * splitCumul[d]);
  }
 
  // Add it to the contained box
  if (index < numBoxes) // avoid segfaults for floating point round-off errors.
    return m_Children[index]->getBoxAtCoord(coords);
  else
    return nullptr;
}
 
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::splitContents(size_t index, Kernel::ThreadScheduler *ts) {
  // You can only split it if it is a MDBox (not MDGridBox).
  MDBox<MDE, nd> *box = dynamic_cast<MDBox<MDE, nd> *>(m_Children[index]);
  if (!box)
    return;
  // Track how many MDBoxes there are in the overall workspace
  this->m_BoxController->trackNumBoxes(box->getDepth());
  // Construct the grid box. This should take the object out of the disk MRU
  auto gridbox = new MDGridBox<MDE, nd>(box);
 
  // Delete the old ungridded box
  delete m_Children[index];
  // And now we have a gridded box instead of a boring old regular box.
  m_Children[index] = gridbox;
 
  if (ts) {
    // Create a task to split the newly created MDGridBox.
    ts->push(std::make_shared<Kernel::FunctionTask>(std::bind(&MDGridBox<MDE, nd>::splitAllIfNeeded, &*gridbox, ts)));
  } else {
    gridbox->splitAllIfNeeded(nullptr);
  }
}
 
//-----------------------------------------------------------------------------------------------
TMDE(size_t MDGridBox)::getChildIndexFromID(size_t childId) const {
  for (size_t index = 0; index < numBoxes; index++) {
    if (m_Children[index]->getID() == childId)
      return index;
  }
  return UNDEF_SIZET;
}
 
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::splitAllIfNeeded(Kernel::ThreadScheduler *ts) {
  for (size_t i = 0; i < numBoxes; ++i) {
    MDBox<MDE, nd> *box = dynamic_cast<MDBox<MDE, nd> *>(m_Children[i]);
    if (box) {
      // Plain MD-Box. Does it need to be split?
      if (this->m_BoxController->willSplit(box->getNPoints(), box->getDepth())) {
        // The MDBox needs to split into a grid box.
        if (!ts) {
          // ------ Perform split serially (no ThreadPool) ------
          auto gridBox = new MDGridBox<MDE, nd>(box);
          // Track how many MDBoxes there are in the overall workspace
          this->m_BoxController->trackNumBoxes(box->getDepth());
          // Replace in the array
          m_Children[i] = gridBox;
          // Delete the old box
          delete box;
          // Now recursively check if this NEW grid box's contents should be
          // split too
          gridBox->splitAllIfNeeded(nullptr);
        } else {
          // ------ Perform split in parallel (using ThreadPool) ------
          // So we create a task to split this MDBox,
          // Task is : this->splitContents(i, ts);
          ts->push(
              std::make_shared<Kernel::FunctionTask>(std::bind(&MDGridBox<MDE, nd>::splitContents, &*this, i, ts)));
        }
      } else {
        // This box does NOT have enough events to be worth splitting, if it do
        // have at least something in memory then,
        Kernel::ISaveable *const pSaver(box->getISaveable());
        if (pSaver && box->getDataInMemorySize() > 0) {
          // Mark the box as "to-write" in DiskBuffer. If the buffer is full,
          // the boxes will be dropped on disk
 
          this->m_BoxController->getFileIO()->toWrite(pSaver);
        }
      }
    } else {
      // It should be a MDGridBox
      MDGridBox<MDE, nd> *gridBox = dynamic_cast<MDGridBox<MDE, nd> *>(m_Children[i]);
      if (gridBox) {
        // Now recursively check if this old grid box's contents should be split
        // too
        if (!ts || (this->nPoints < this->m_BoxController->getAddingEvents_eventsPerTask()))
          // Go serially if there are only a few points contained (less
          // overhead).
          gridBox->splitAllIfNeeded(ts);
        else
          // Go parallel if this is a big enough gridbox.
          // Task is : gridBox->splitAllIfNeeded(ts);
          ts->push(
              std::make_shared<Kernel::FunctionTask>(std::bind(&MDGridBox<MDE, nd>::splitAllIfNeeded, &*gridBox, ts)));
      }
    }
  }
}
 
TMDE(void MDGridBox)::integrateSphere(API::CoordTransform &radiusTransform, const coord_t radiusSquared,
                                      signal_t &signal, signal_t &errorSquared, const coord_t innerRadiusSquared,
                                      const bool useOnePercentBackgroundCorrection) const {
  // We start by looking at the vertices at every corner of every box contained,
  // to see which boxes are partially contained/fully contained.
 
  // One entry with the # of vertices in this box contained; start at 0.
  std::vector<size_t> verticesContained(numBoxes, 0);
 
  // Set to true if there is a possibility of the box at least partly touching
  // the integration volume.
  std::vector<bool> boxMightTouch(numBoxes, 0);
 
  // How many vertices does one box have? 2^nd, or bitwise shift left 1 by nd
  // bits
  size_t maxVertices = 1 << nd;
 
  // set up caches for box sizes and min box values
  coord_t boxSize[nd];
  coord_t minBoxVal[nd];
 
  // The number of vertices in each dimension is the # split[d] + 1
  size_t vertices_max[nd];
  Kernel::Utils::NestedForLoop::SetUp(nd, vertices_max, 0);
  for (size_t d = 0; d < nd; ++d) {
    vertices_max[d] = split[d] + 1;
    // cache box sizes and min box valyes for performance
    boxSize[d] = static_cast<coord_t>(m_SubBoxSize[d]);
    minBoxVal[d] = static_cast<coord_t>(this->extents[d].getMin());
  }
 
  // The index to the vertex in each dimension
  size_t vertexIndex[nd];
  Kernel::Utils::NestedForLoop::SetUp(nd, vertexIndex, 0);
  size_t boxIndex[nd];
  Kernel::Utils::NestedForLoop::SetUp(nd, boxIndex, 0);
  size_t indexMaker[nd];
  Kernel::Utils::NestedForLoop::SetUpIndexMaker(nd, indexMaker, split);
 
  bool allDone = false;
  while (!allDone) {
    // Coordinates of this vertex
    coord_t vertexCoord[nd];
    for (size_t d = 0; d < nd; ++d)
      vertexCoord[d] = static_cast<coord_t>(vertexIndex[d]) * boxSize[d] + minBoxVal[d];
    // static_cast<coord_t>(vertexIndex[d]) * boxSize[d] + this->extents[d].min
 
    // Is this vertex contained?
    coord_t out[nd];
    radiusTransform.apply(vertexCoord, out);
    if (out[0] < radiusSquared && out[0] > innerRadiusSquared) {
      // Yes, this vertex is contained within the integration volume!
      //        std::cout << "vertex at " << vertexCoord[0] << ", " <<
      //        vertexCoord[1] << ", " << vertexCoord[2] << " is contained\n";
 
      // This vertex is shared by up to 2^nd adjacent boxes (left-right along
      // each dimension).
      for (size_t neighb = 0; neighb < maxVertices; ++neighb) {
        // The index of the box is the same as the vertex, but maybe - 1 in each
        // possible combination of dimensions
        bool badIndex = false;
        // Build the index of the neighbor
        for (size_t d = 0; d < nd; d++) {
          boxIndex[d] = vertexIndex[d] - ((neighb & ((size_t)1 << d)) >> d); //(this does a bitwise and mask,
          // shifted back to 1 to subtract 1
          // to the dimension)
          // Taking advantage of the fact that unsigned(0)-1 = some large
          // POSITIVE number.
          if (boxIndex[d] >= split[d]) {
            badIndex = true;
            break;
          }
        }
        if (!badIndex) {
          // Convert to linear index
          size_t linearIndex = Kernel::Utils::NestedForLoop::GetLinearIndex(nd, boxIndex, indexMaker);
          // So we have one more vertex touching this box that is contained in
          // the integration volume. Whew!
          verticesContained[linearIndex]++;
          //            std::cout << "... added 1 vertex to box " <<
          //            boxes[linearIndex]->getExtentsStr() << "\n";
        }
      }
    }
 
    // Increment the counter(s) in the nested for loops.
    allDone = Kernel::Utils::NestedForLoop::Increment(nd, vertexIndex, vertices_max);
  }
 
  // NOTE:
  //    The following section is just trying to uncover the peak center so that
  //    we can compute the distance bewteen peak center and box center correctly
  //    without worrying about the skew coming from coordTransformDistnace.
  auto tmpRadiusTransform = dynamic_cast<CoordTransformDistance *>(&radiusTransform);
  if (tmpRadiusTransform == nullptr) {
    throw std::runtime_error("radiusTransform has to be CoordTransformDistance");
  }
  auto peakCenter = tmpRadiusTransform->getCenter();
  double peakRadius = std::sqrt(radiusSquared);
  double peakInnerRadius = std::sqrt(innerRadiusSquared);
 
  // OK, we've done counting all the vertices.
  // Now let's go through and check each box.
  for (size_t bIndex = 0; bIndex < numBoxes; ++bIndex) {
    API::IMDNode *box = m_Children[bIndex];
 
    // First, check if we have reached the base case where the box is completely
    // enveloped by the peak
    if (verticesContained[bIndex] >= maxVertices) {
      // Use the integrated sum of signal in the box
      signal += box->getSignal();
      errorSquared += box->getErrorSquared();
      continue; // move on to next
    }
 
    // Second, check if there is at least one vertex in the integration volume,
    // and kick off recursive search until we reach the base case
    if (verticesContained[bIndex] > 0) {
      box->integrateSphere(radiusTransform, radiusSquared, signal, errorSquared, innerRadiusSquared,
                           useOnePercentBackgroundCorrection);
      continue;
    }
 
    // Last, no vertices in the integration volume
    //  -- there are only two cases (see below) where we can
    //     skip the box and its children, and it requires
    //     the knowledge of how box and peak are spatially
    //     positioned in Euclidian space.
    // NOTE:
    //  For ellipsoid, we are using the long semi-axis, which
    //  means we will inevitablly search some empty boxes, but
    //  this is by design as we do not want to miss any box
    //  that might contain events we need to collect.
    coord_t boxCenter[nd];
    box->getCenter(boxCenter);
    double distPeakCenterToBoxCenter = 0.0;
    for (size_t i = 0; i < nd; ++i) {
      distPeakCenterToBoxCenter += (boxCenter[i] - peakCenter[i]) * (boxCenter[i] - peakCenter[i]);
    }
    distPeakCenterToBoxCenter = std::sqrt(distPeakCenterToBoxCenter);
    double boxRadius = std::sqrt(diagonalSquared);
 
    //  - Box is completely isolated from peak
    //        distPeakCenterToBoxCenter - peakRadius > boxRadius
    //    -> stop search and move on
    if (distPeakCenterToBoxCenter - peakRadius > boxRadius) {
      // Debug output
      // std::ostringstream debugmsg;
      // debugmsg << "R_peak = " << peakRadius << "\n"
      //          << "R_box = " << boxRadius << "\n"
      //          << "Peak Center: ";
      // for (size_t i = 0; i < nd; ++i) {
      //   debugmsg << peakCenter[i] << ",";
      // }
      // debugmsg << "\n"
      //          << "Box center: ";
      // for (size_t i = 0; i < nd; ++i) {
      //   debugmsg << boxCenter[i] << ",";
      // }
      // debugmsg << "\n"
      //          << "distPeakCenterToBoxCenter = " << distPeakCenterToBoxCenter
      //          << "\n";
      // std::cout << debugmsg.str();
      continue;
    }
    //  - Tiny box falls into the donut hole
    //        distPeakCenterToBoxCenter + boxRadius < peakInnerRadius
    //    -> stop search and move on
    if (peakInnerRadius > 0 && distPeakCenterToBoxCenter + boxRadius < peakInnerRadius) {
      continue;
    }
    //  - All other cases, we need to refine the box, i.e. recursively
    //    search the child box
    box->integrateSphere(radiusTransform, radiusSquared, signal, errorSquared, innerRadiusSquared,
                         useOnePercentBackgroundCorrection);
  } // (for each box)
}
 
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::centroidSphere(API::CoordTransform &radiusTransform, const coord_t radiusSquared,
                                     coord_t *centroid, signal_t &signal) const {
  for (size_t i = 0; i < numBoxes; ++i) {
    // Go through each contained box
    API::IMDNode *box = m_Children[i];
    coord_t boxCenter[nd];
    box->getCenter(boxCenter);
 
    // Distance from center to the peak integration center
    coord_t out[nd];
    radiusTransform.apply(boxCenter, out);
 
    if (out[0] < diagonalSquared * 0.72 + radiusSquared) {
      // If the center is closer than the size of the box, then it MIGHT be
      // touching.
      // (We multiply by 0.72 (about sqrt(2)) to look for half the diagonal).
      // NOTE! Watch out for non-spherical transforms!
 
      // Go down one level to keep centroiding
      box->centroidSphere(radiusTransform, radiusSquared, centroid, signal);
    }
  } // (for each box)
}
//-----------------------------------------------------------------------------------------------
GNU_DIAG_OFF("array-bounds")
TMDE(void MDGridBox)::integrateCylinder(Mantid::API::CoordTransform &radiusTransform, const coord_t radius,
                                        const coord_t length, signal_t &signal, signal_t &errorSquared,
                                        std::vector<signal_t> &signal_fit) const {
  // We start by looking at the vertices at every corner of every box contained,
  // to see which boxes are partially contained/fully contained.
 
  // One entry with the # of vertices in this box contained; start at 0.
  auto verticesContained = new size_t[numBoxes];
  memset(verticesContained, 0, numBoxes * sizeof(size_t));
 
  // Set to true if there is a possibility of the box at least partly touching
  // the integration volume.
  auto boxMightTouch = new bool[numBoxes];
  memset(boxMightTouch, 0, numBoxes * sizeof(bool));
 
  // How many vertices does one box have? 2^nd, or bitwise shift left 1 by nd
  // bits
  size_t maxVertices = 1 << nd;
 
  // set up caches for box sizes and min box values
  coord_t boxSize[nd];
  coord_t minBoxVal[nd];
 
  // The number of vertices in each dimension is the # split[d] + 1
  size_t vertices_max[nd];
  Kernel::Utils::NestedForLoop::SetUp(nd, vertices_max, 0);
  for (size_t d = 0; d < nd; ++d) {
    vertices_max[d] = split[d] + 1;
    // cache box sizes and min box valyes for performance
    boxSize[d] = static_cast<coord_t>(m_SubBoxSize[d]);
    minBoxVal[d] = static_cast<coord_t>(this->extents[d].getMin());
  }
 
  // The index to the vertex in each dimension
  size_t vertexIndex[nd];
  Kernel::Utils::NestedForLoop::SetUp(nd, vertexIndex, 0);
  size_t boxIndex[nd];
  Kernel::Utils::NestedForLoop::SetUp(nd, boxIndex, 0);
  size_t indexMaker[nd];
  Kernel::Utils::NestedForLoop::SetUpIndexMaker(nd, indexMaker, split);
 
  size_t numSteps = signal_fit.size();
  double deltaQ = length / static_cast<double>(numSteps - 1);
  bool allDone = false;
  while (!allDone) {
    // Coordinates of this vertex
    coord_t vertexCoord[nd];
    for (size_t d = 0; d < nd; ++d)
      vertexCoord[d] = static_cast<coord_t>(vertexIndex[d]) * boxSize[d] + minBoxVal[d];
    // static_cast<coord_t>(vertexIndex[d]) * boxSize[d] + this->extents[d].min
 
    // Is this vertex contained?
    coord_t out[2]; // radius and length of cylinder
    radiusTransform.apply(vertexCoord, out);
    if (out[0] < radius && std::fabs(out[1]) < 0.5 * length) {
      // Yes, this vertex is contained within the integration volume!
      //        std::cout << "vertex at " << vertexCoord[0] << ", " <<
      //        vertexCoord[1] << ", " << vertexCoord[2] << " is contained\n";
 
      // This vertex is shared by up to 2^nd adjacent boxes (left-right along
      // each dimension).
      for (size_t neighb = 0; neighb < maxVertices; ++neighb) {
        // The index of the box is the same as the vertex, but maybe - 1 in each
        // possible combination of dimensions
        bool badIndex = false;
        // Build the index of the neighbor
        for (size_t d = 0; d < nd; d++) {
          boxIndex[d] = vertexIndex[d] - ((neighb & ((size_t)1 << d)) >> d); //(this does a bitwise and mask,
          // shifted back to 1 to subtract 1
          // to the dimension)
          // Taking advantage of the fact that unsigned(0)-1 = some large
          // POSITIVE number.
          if (boxIndex[d] >= split[d]) {
            badIndex = true;
            break;
          }
        }
        if (!badIndex) {
          // Convert to linear index
          size_t linearIndex = Kernel::Utils::NestedForLoop::GetLinearIndex(nd, boxIndex, indexMaker);
          // So we have one more vertex touching this box that is contained in
          // the integration volume. Whew!
          verticesContained[linearIndex]++;
          //            std::cout << "... added 1 vertex to box " <<
          //            boxes[linearIndex]->getExtentsStr() << "\n";
        }
      }
    }
 
    // Increment the counter(s) in the nested for loops.
    allDone = Kernel::Utils::NestedForLoop::Increment(nd, vertexIndex, vertices_max);
  }
 
  // OK, we've done all the vertices. Now we go through and check each box.
  for (size_t i = 0; i < numBoxes; ++i) {
    API::IMDNode *box = m_Children[i];
    // Box partially contained?
    bool partialBox = false;
 
    // Is this box fully contained?
    if (verticesContained[i] >= maxVertices) {
      std::vector<coord_t> coordTable;
      size_t nColumns;
      box->getEventsData(coordTable, nColumns);
      if (nColumns > 0 && nd > 1) {
        size_t nEvents = coordTable.size() / nColumns;
        size_t skipCol = 2; // lean events
        if (nColumns == 7)
          skipCol += 2; // events
        for (size_t k = 0; k < nEvents; k++) {
          coord_t eventCenter[nd];
          for (size_t l = 0; l < nd; l++)
            eventCenter[l] = coordTable[k * nColumns + skipCol + l];
          coord_t out[nd];
          radiusTransform.apply(eventCenter, out);
          // add event to appropriate y channel
          size_t xchannel = static_cast<size_t>(std::floor(out[1] / deltaQ)) + numSteps / 2;
 
          if (xchannel < numSteps)
            signal_fit[xchannel] += coordTable[k * nColumns];
        }
      }
      // box->releaseEvents();
      // Use the integrated sum of signal in the box
      signal += box->getSignal();
      errorSquared += box->getErrorSquared();
 
      //        std::cout << "box at " << i << " (" << box->getExtentsStr() <<
      //        ") is fully contained. Vertices = " << verticesContained[i] <<
      //        "\n";
      // Go on to the next box
      continue;
    }
 
    if (verticesContained[i] == 0) {
      // There is a chance that this part of the box is within integration
      // volume,
      // even if no vertex of it is.
      coord_t boxCenter[nd];
      box->getCenter(boxCenter);
 
      // Distance from center to the peak integration center
      coord_t out[nd];
      radiusTransform.apply(boxCenter, out);
      if ((nd >= 1) && out[0] < std::sqrt(diagonalSquared * 0.72 + radius * radius) &&
          (nd >= 2 && std::fabs(out[1]) < std::sqrt(diagonalSquared * 0.72 + 0.25 * length * length))) {
        // If the center is closer than the size of the box, then it MIGHT be
        // touching.
        // (We multiply by 0.72 (about sqrt(2)) to look for half the diagonal).
        // NOTE! Watch out for non-spherical transforms!
        //          std::cout << "box at " << i << " is maybe touching\n";
        partialBox = true;
      }
    } else {
      partialBox = true;
      //        std::cout << "box at " << i << " has a vertex touching\n";
    }
 
    // We couldn't rule out that the box might be partially contained.
    if (partialBox) {
      // Use the detailed integration method.
      box->integrateCylinder(radiusTransform, radius, length, signal, errorSquared, signal_fit);
      //        std::cout << ".signal=" << signal << "\n";
    }
  } // (for each box)
 
  //    std::cout << "Depth " << this->getDepth() << " with " <<
  //    numFullyContained << " fully contained; " << numPartiallyContained << "
  //    partial. Signal = " << signal <<"\n";
 
  delete[] verticesContained;
  delete[] boxMightTouch;
}
GNU_DIAG_ON("array-bounds")
 
 
TMDE(bool MDGridBox)::getIsMasked() const {
  bool isMasked = false;
  for (size_t i = 0; i < numBoxes; ++i) {
    // Go through each contained box
    API::IMDNode *box = m_Children[i];
    if (box->getIsMasked()) {
      isMasked = true;
      break;
    }
  }
  return isMasked;
}
 
TMDE(void MDGridBox)::mask() {
  for (size_t i = 0; i < numBoxes; ++i) {
    // Go through each contained box
    API::IMDNode *box = m_Children[i];
    box->mask();
  }
}
 
TMDE(void MDGridBox)::unmask() {
  for (size_t i = 0; i < numBoxes; ++i) {
    // Go through each contained box
    API::IMDNode *box = m_Children[i];
    box->unmask();
  }
}
//------------------------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------------------------
/* Internal TMP class to simplify adding events to the box for events and lean
 * events using single interface. One would nead to overload the box class
 * otherwise*/
template <typename MDE, size_t nd> struct IF_EVENT {
public:
  // create generic events from array of events data and add them to the grid
  // box
  static inline void EXEC(MDGridBox<MDE, nd> *pBox, const std::vector<signal_t> &sigErrSq,
                          const std::vector<coord_t> &Coord, const std::vector<uint16_t> &expInfoIndex,
                          const std::vector<uint16_t> &goniometerIndex, const std::vector<uint32_t> &detectorId,
                          size_t nEvents) {
    for (size_t i = 0; i < nEvents; i++)
      pBox->addEvent(MDEvent<nd>(sigErrSq[2 * i], sigErrSq[2 * i + 1], expInfoIndex[i], goniometerIndex[i],
                                 detectorId[i], &Coord[i * nd]));
  }
};
/* Specialize for the case of LeanEvent */
template <size_t nd> struct IF_EVENT<MDLeanEvent<nd>, nd> {
public:
  // create lean events from array of events data and add them to the grid box
  static inline void EXEC(MDGridBox<MDLeanEvent<nd>, nd> *pBox, const std::vector<signal_t> &sigErrSq,
                          const std::vector<coord_t> &Coord, const std::vector<uint16_t> & /*expInfoIndex*/,
                          const std::vector<uint16_t> & /*goniometerIndex*/,
                          const std::vector<uint32_t> & /*detectorId*/, size_t nEvents) {
    for (size_t i = 0; i < nEvents; i++)
      pBox->addEvent(MDLeanEvent<nd>(sigErrSq[2 * i], sigErrSq[2 * i + 1], &Coord[i * nd]));
  }
};
 
TMDE(size_t MDGridBox)::buildAndAddEvents(const std::vector<signal_t> &sigErrSq, const std::vector<coord_t> &Coord,
                                          const std::vector<uint16_t> &expInfoIndex,
                                          const std::vector<uint16_t> &goniometerIndex,
                                          const std::vector<uint32_t> &detectorId) {
 
  size_t nEvents = sigErrSq.size() / 2;
  IF_EVENT<MDE, nd>::EXEC(this, sigErrSq, Coord, expInfoIndex, goniometerIndex, detectorId, nEvents);
 
  return 0;
}
 
TMDE(void MDGridBox)::buildAndAddEvent(const signal_t Signal, const signal_t errorSq, const std::vector<coord_t> &point,
                                       uint16_t expInfoIndex, uint16_t goniometerIndex, uint32_t detectorId) {
  this->addEvent(IF<MDE, nd>::BUILD_EVENT(Signal, errorSq, &point[0], expInfoIndex, goniometerIndex, detectorId));
}
 
//-----------------------------------------------------------------------------------------------
TMDE(void MDGridBox)::buildAndAddEventUnsafe(const signal_t Signal, const signal_t errorSq,
                                             const std::vector<coord_t> &point, uint16_t expInfoIndex,
                                             uint16_t goniometerIndex, uint32_t detectorId) {
  this->addEventUnsafe(IF<MDE, nd>::BUILD_EVENT(Signal, errorSq, &point[0], expInfoIndex, goniometerIndex, detectorId));
}
 
//-----------------------------------------------------------------------------------------------
TMDE(inline size_t MDGridBox)::addEvent(const MDE &event) {
  size_t cindex = calculateChildIndex(event);
 
  // We can erroneously get cindex == numBoxes for events which fall on the
  // upper boundary of the last child box, so add these events to the last box
  if (cindex == numBoxes)
    cindex = numBoxes - 1;
 
  if (cindex < numBoxes)
    return m_Children[cindex]->addEvent(event);
  else
    return 0;
}
 
//-----------------------------------------------------------------------------------------------
TMDE(inline size_t MDGridBox)::addEventUnsafe(const MDE &event) {
  size_t cindex = calculateChildIndex(event);
 
  // We can erroneously get cindex == numBoxes for events which fall on the
  // upper boundary of the last child box, so add these events to the last box
  if (cindex == numBoxes)
    cindex = numBoxes - 1;
 
  if (cindex < numBoxes)
    return m_Children[cindex]->addEventUnsafe(event);
  else
    return 0;
}
 
TMDE(inline void MDGridBox)::setChild(size_t index, MDGridBox<MDE, nd> *newChild) {
  // Delete the old box  (supposetly ungridded);
  delete this->m_Children[index];
  // set new box, supposetly gridded
  this->m_Children[index] = newChild;
}
TMDE(void MDGridBox)::setFileBacked(const uint64_t /*fileLocation*/, const size_t /*fileSize*/,
                                    const bool /*markSaved*/) {
  throw(Kernel::Exception::NotImplementedError("Recursive file backed is not "
                                               "yet implemented (unclear how "
                                               "to set file location etc)"));
}
TMDE(void MDGridBox)::setFileBacked() {
  for (size_t i = 0; i < this->numBoxes; i++) {
    m_Children[i]->setFileBacked();
  }
}
TMDE(void MDGridBox)::clearFileBacked(bool loadDiskBackedData) {
  auto it = m_Children.begin();
  auto it_end = m_Children.end();
  for (; it != it_end; it++) {
    (*it)->clearFileBacked(loadDiskBackedData);
  }
}
 
TMDE(size_t MDGridBox)::calculateChildIndex(const MDE &event) const {
  size_t cindex(0);
  for (size_t d = 0; d < nd; d++) {
    // Accumulate the index
    const auto coordinate = event.getCenter(d);
    const auto offset = coordinate - this->extents[d].getMin();
    const int splitD = static_cast<int>(split[d]);
    // clamp to splitD-1: FP rounding can place an in-box coordinate one slot past the last valid child
    auto childIndex = static_cast<int>(offset * m_SubBoxSizeInv[d]);
    if (childIndex == splitD && coordinate <= this->extents[d].getMax())
      childIndex = splitD - 1;
    cindex += childIndex * splitCumul[d];
  }
  return cindex;
}
} // namespace DataObjects
 
} // namespace Mantid