MDGridBox.tcc

#include "MantidKernel/Task.h"
#include "MantidKernel/Utils.h"
#include "MantidKernel/FunctionTask.h"
#include "MantidKernel/Timer.h"
#include "MantidKernel/ThreadPool.h"
#include "MantidKernel/ThreadScheduler.h"
#include "MantidKernel/ThreadSchedulerMutexes.h"
#include "MantidKernel/WarningSuppressions.h"
#include "MantidDataObjects/MDBox.h"
#include "MantidDataObjects/MDEvent.h"
#include "MantidDataObjects/MDGridBox.h"
#include <boost/math/special_functions/round.hpp>
#include <boost/optional.hpp>
#include <ostream>
#include "MantidKernel/Strings.h"

// 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 {

////===============================================================================================
////===============================================================================================

//-----------------------------------------------------------------------------------------------
/** Constructor with a box controller.
 * @param bc :: poineter to the BoxController, owned by workspace
 * @param depth :: recursive split depth
 * @param extentsVector :: size of the box
 */
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();
}

/** convenience Constructor, taking the shared pointer and extracting const
 * pointer from it
  * @param bc :: shared poineter to the BoxController, owned by workspace
  * @param depth :: recursive split depth
  * @param extentsVector :: size of the box
*/
TMDE(MDGridBox)::MDGridBox(
    boost::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();
}
/// common part of MDGridBox contstructor;
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.
  boost::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.get()[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;
}

//-----------------------------------------------------------------------------------------------
/** Constructor
 * @param box :: MDBox containing the events to split
 */
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();
}
/**Internal function to do main job of filling in a GridBox contents  (part of
 * the constructor) */
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.push_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
}

//-----------------------------------------------------------------------------------------------
/** Copy constructor
 * @param other :: MDGridBox to copy
 * @param otherBC :: mandatory pointer to other box controller, which will split
 this box.
                     if it the same BC, as the one for the copied box, it needs
 to be taken explicitly from the
                     copied 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];
  }
  // 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.push_back(newBox);
    } else if (otherMDGridBox) {
      auto newBox = new MDGridBox<MDE, nd>(*otherMDGridBox, otherBC);
      newBox->setParent(this);
      m_Children.push_back(newBox);
    } else {
      throw std::runtime_error(
          "MDGridBox::copy_ctor(): an unexpected child box type was found.");
    }
  }
}

//-----------------------------------------------------------------------------------------------
/** Transform the dimensions contained in this box
 * x' = x*scaling + offset.
 * NON-RECURSIVE!
 *
 * @param scaling :: multiply each coordinate by this value.
 * @param offset :: after multiplying, add this offset.
 */
TMDE(void MDGridBox)::transformDimensions(std::vector<double> &scaling,
                                          std::vector<double> &offset) {
  MDBoxBase<MDE, nd>::transformDimensions(scaling, offset);
  this->computeSizesFromSplit();
}

//-----------------------------------------------------------------------------------------------
/** Compute some data from the split[] array and the extents.
 *
 * @return :: the total number of boxes */
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]);
    // Accumulate the squared diagonal length.
    diagSum += m_SubBoxSize[d] * m_SubBoxSize[d];
  }
  diagonalSquared = static_cast<coord_t>(diagSum);

  return tot;
}

//-----------------------------------------------------------------------------------------------
/// Destructor
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();
}

//-----------------------------------------------------------------------------------------------
/** Clear any points contained. */
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();
  }
}

//-----------------------------------------------------------------------------------------------
/** Returns the number of dimensions in this box */
TMDE(size_t MDGridBox)::getNumDims() const { return nd; }

//-----------------------------------------------------------------------------------------------
/// Recursiveluy calculates the amount of the data located in memory. Slow
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;
}

//-----------------------------------------------------------------------------------------------
/** Returns the number of un-split MDBoxes in this box (recursively including
 * all children)
 * @return :: the total # of MDBoxes in all children */
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;
}

//-----------------------------------------------------------------------------------------------
/** @return The number of children of the MDGridBox, not recursively
 */
TMDE(size_t MDGridBox)::getNumChildren() const { return numBoxes; }

//-----------------------------------------------------------------------------------------------
/** Get a child box
 * @param index :: index into the array, within range 0..getNumChildren()-1
 * @return the child MDBoxBase pointer.
 */
TMDE(API::IMDNode *MDGridBox)::getChild(size_t index) {
  return m_Children[index];
}

//-----------------------------------------------------------------------------------------------
/** Directly set the children of the MDGridBox. Used in file loading.
 * Should not be called on a box with children; the existing children are NOT
 *deleted.
 *
 * @param otherBoxes:: reference to a vector of boxes containing the children
 * @param indexStart :: start point in the vector
 * @param indexEnd :: end point in the vector, not-inclusive
 */
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.push_back(dynamic_cast<MDBoxBase<MDE, nd> *>(*it));
    m_Children.back()->setParent(this);
  }
  numBoxes = m_Children.size();
}

//-----------------------------------------------------------------------------------------------
/** Helper function to get the index into the linear array given
 * an array of indices for each dimension (0 to nd)
 * @param indices :: array of size[nd]
 * @return size_t index into m_Children[].
 */
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;
}

//-----------------------------------------------------------------------------------------------
/** Refresh the cache of nPoints, signal and error,
 * by adding up all boxes (recursively).
 * MDBoxes' totals are used directly.
 *
 * @param ts :: ThreadScheduler pointer to perform the caching
 *  in parallel. If NULL, it will be performed in series.
 */
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");
  }
}

//-----------------------------------------------------------------------------------------------
/** Allocate and return a vector with a copy of all events contained
 */
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;
}

//-----------------------------------------------------------------------------------------------
/** Return all boxes contained within.
 *
 * @param outBoxes :: vector to fill
 * @param maxDepth :: max depth value of the returned boxes.
 * @param leafOnly :: if true, only add the boxes that are no more subdivided
 *(leaves on the tree)
 */
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.push_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.push_back(this);
  }
}

//-----------------------------------------------------------------------------------------------
/** Return all boxes contained within, limited by an implicit function.
 *
 * This method evaluates each vertex to see how it is contained by the implicit
 *function.
 * For example, if there are 4x4 boxes, there are 5x5 vertices to evaluate.
 *
 * All boxes that might be touching the implicit function are returned
 *(including ones that
 * overlap without any point actually in the function).
 *
 * @param outBoxes :: vector to fill
 * @param maxDepth :: max depth value of the returned boxes.
 * @param leafOnly :: if true, only add the boxes that are no more subdivided
 *(leaves on the tree)
 * @param function :: implicitFunction pointer
 */
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.push_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.push_back(this);
  }
}

//-----------------------------------------------------------------------------------------------
/** Returns the lowest-level box at the given coordinates
 * @param coords :: nd-sized array of the coordinate of the point to look at
 * @return MDBoxBase pointer.
 */
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;
}

//-----------------------------------------------------------------------------------------------
/** Split a box that is contained in the GridBox, at the given index,
 * into a MDGridBox.
 *
 * Thread-safe as long as 'index' is different for all threads.
 *
 * @param index :: index into the boxes vector.
 *        Warning: No bounds check is made, don't give stupid values!
 * @param ts :: optional ThreadScheduler * that will be used to parallelize
 *        recursive splitting. Set to NULL for no recursive splitting.
 */
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(new Kernel::FunctionTask(
        boost::bind(&MDGridBox<MDE, nd>::splitAllIfNeeded, &*gridbox, ts)));
  } else {
    gridbox->splitAllIfNeeded(nullptr);
  }
}

//-----------------------------------------------------------------------------------------------
/** Get the child index from its ID
 *
 * @param childId :: ID of the child you want
 * @return the index into the children of this grid box; size_t(-1) if NOT
 *found.
 */
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;
}

//-----------------------------------------------------------------------------------------------
/** Goes through all the sub-boxes and splits them if they contain
 * enough events to be worth it.
 *
 * @param ts :: optional ThreadScheduler * that will be used to parallelize
 *        recursive splitting. Set to NULL to do it serially.
 */
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(new Kernel::FunctionTask(
              boost::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(new Kernel::FunctionTask(boost::bind(
              &MDGridBox<MDE, nd>::splitAllIfNeeded, &*gridBox, ts)));
      }
    }
  }
}

//-----------------------------------------------------------------------------------------------
/** Perform centerpoint binning of events, with bins defined
 * in axes perpendicular to the axes of the workspace.
 *
 * @param bin :: MDBin object giving the limits of events to accept.
 * @param fullyContained :: optional bool array sized [nd] of which dimensions
 *are known to be fully contained (for MDSplitBox)
 */
TMDE(void MDGridBox)::centerpointBin(MDBin<MDE, nd> &bin,
                                     bool *fullyContained) const {

  // The MDBin ranges from index_min to index_max (inclusively) if each
  // dimension. So
  // we'll need to make nested loops from index_min[0] to index_max[0]; from
  // index_min[1] to index_max[1]; etc.
  int index_min[nd];
  int index_max[nd];
  // For running the nested loop, counters of each dimension. These are bounded
  // by 0..split[d]
  size_t counters_min[nd];
  size_t counters_max[nd];

  for (size_t d = 0; d < nd; d++) {
    int min, max;

    // The min index in this dimension (we round down - we'll include this edge)
    if (bin.m_min[d] >= this->extents[d].getMin()) {
      min = int((bin.m_min[d] - this->extents[d].getMin()) / m_SubBoxSize[d]);
      counters_min[d] = min;
    } else {
      min = -1; // Goes past the edge
      counters_min[d] = 0;
    }

    // If the minimum is bigger than the number of blocks in that dimension,
    // then the bin is off completely in
    //  that dimension. There is nothing to integrate.
    if (min >= static_cast<int>(split[d]))
      return;
    index_min[d] = min;

    // The max index in this dimension (we round UP, but when we iterate we'll
    // NOT include this edge)
    if (bin.m_max[d] < this->extents[d].getMax()) {
      max = int(ceil((bin.m_max[d] - this->extents[d].getMin()) /
                     m_SubBoxSize[d])) -
            1;
      counters_max[d] =
          max + 1; // (the counter looping will NOT include counters_max[d])
    } else {
      max = int(split[d]);   // Goes past THAT edge
      counters_max[d] = max; // (the counter looping will NOT include max)
    }

    // If the max value is before the min, that means NOTHING is in the bin, and
    // we can return
    if ((max < min) || (max < 0))
      return;
    index_max[d] = max;

    // std::cout << d << " from " << std::setw(5) << index_min[d] << " to " <<
    // std::setw(5)  << index_max[d] << "inc\n";
  }

  // If you reach here, than at least some of bin is overlapping this box
  size_t counters[nd];
  for (size_t d = 0; d < nd; d++)
    counters[d] = counters_min[d];

  bool allDone = false;
  while (!allDone) {
    size_t index = getLinearIndex(counters);
    // std::cout << index << ": " << counters[0] << ", " << counters[1] <<
    // '\n';

    // Find if the box is COMPLETELY held in the bin.
    bool completelyWithin = true;
    for (size_t dim = 0; dim < nd; dim++)
      if ((static_cast<int>(counters[dim]) <= index_min[dim]) ||
          (static_cast<int>(counters[dim]) >= index_max[dim])) {
        // The index we are at is at the edge of the integrated area (index_min
        // or index_max-1)
        // That means that the bin only PARTIALLY covers this MDBox
        completelyWithin = false;
        break;
      }

    if (completelyWithin) {
      // Box is completely in the bin.
      // std::cout << "Box at index " << counters[0] << ", " << counters[1] << "
      // is entirely contained.\n";
      // Use the aggregated signal and error
      bin.m_signal += m_Children[index]->getSignal();
      bin.m_errorSquared += m_Children[index]->getErrorSquared();
    } else {
      // Perform the binning
      m_Children[index]->centerpointBin(bin, fullyContained);
    }

    // Increment the counter(s) in the nested for loops.
    allDone = Kernel::Utils::NestedForLoop::Increment(
        nd, counters, counters_max, counters_min);
  }
}

//  TMDE(
//  void MDGridBox)::generalBin(MDBin<MDE,nd> & bin,
//  Mantid::API::ImplicitFunction & function) const
//  {
//    // The MDBin ranges from index_min to index_max (inclusively) if each
//    dimension. So
//    // we'll need to make nested loops from index_min[0] to index_max[0]; from
//    index_min[1] to index_max[1]; etc.
//    int index_min[nd];
//    int index_max[nd];
//    // For running the nested loop, counters of each dimension. These are
//    bounded by 0..split[d]
//    size_t counters_min[nd];
//    size_t counters_max[nd];
//
//    for (size_t d=0; d<nd; d++)
//    {
//      int min,max;
//
//      // The min index in this dimension (we round down - we'll include this
//      edge)
//      if (bin.m_min[d] >= this->extents[d].getMin())
//      {
//        min = int((bin.m_min[d] - this->extents[d].getMin()) / boxSize[d]);
//        counters_min[d] = min;
//      }
//      else
//      {
//        min = -1; // Goes past the edge
//        counters_min[d] = 0;
//      }
//
//      // If the minimum is bigger than the number of blocks in that dimension,
//      then the bin is off completely in
//      //  that dimension. There is nothing to integrate.
//      if (min >= static_cast<int>(split[d]))
//        return;
//      index_min[d] = min;
//
//      // The max index in this dimension (we round UP, but when we iterate
//      we'll NOT include this edge)
//      if (bin.m_max[d] < this->extents[d].max)
//      {
//        max = int(ceil((bin.m_max[d] - this->extents[d].getMin()) /
//        boxSize[d])) - 1;
//        counters_max[d] = max+1; // (the counter looping will NOT include
//        counters_max[d])
//      }
//      else
//      {
//        max = int(split[d]); // Goes past THAT edge
//        counters_max[d] = max; // (the counter looping will NOT include max)
//      }
//
//      // If the max value is before the min, that means NOTHING is in the bin,
//      and we can return
//      if ((max < min) || (max < 0))
//        return;
//      index_max[d] = max;
//
//      //std::cout << d << " from " << std::setw(5) << index_min[d] << " to "
//      << std::setw(5)  << index_max[d] << "inc\n";
//    }
//
//    // If you reach here, than at least some of bin is overlapping this box
//
//
//    // 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.
//    size_t * 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.
//    bool * 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;
//
//    // The index to the vertex in each dimension
//    size_t * vertexIndex = Utils::NestedForLoop::SetUp(nd, 0);
//
//    // This is the index in each dimension at which we start looking at
//    vertices
//    size_t * vertices_min = Utils::NestedForLoop::SetUp(nd, 0);
//    for (size_t d=0; d<nd; ++d)
//    {
//      vertices_min[d] = counters_min[d];
//      vertexIndex[d] = vertices_min[d]; // This is where we start
//    }
//
//    // There is one more vertex in each dimension than there are boxes we are
//    considering
//    size_t * vertices_max = Utils::NestedForLoop::SetUp(nd, 0);
//    for (size_t d=0; d<nd; ++d)
//      vertices_max[d] = counters_max[d]+1;
//
//    size_t * boxIndex = Utils::NestedForLoop::SetUp(nd, 0);
//    size_t * indexMaker = Utils::NestedForLoop::SetUpIndexMaker(nd, split);
//
//    bool allDone = false;
//    while (!allDone)
//    {
//      // Coordinates of this vertex
//      coord_t vertexCoord[nd];
//      bool masks[nd];
//      for (size_t d=0; d<nd; ++d)
//      {
//        vertexCoord[d] = double(vertexIndex[d]) * boxSize[d] +
//        this->extents[d].getMin();
//        masks[d] = false; //HACK ... assumes that all vertexes are used.
//      }
//      // Is this vertex contained?
//      if (function.evaluate(vertexCoord, masks, nd))
//      {
//        // 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 & (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;