MDGridBoxTest.h

#ifndef MDGRIDBOXTEST_H
#define MDGRIDBOXTEST_H

#include "MantidGeometry/MDGeometry/MDBoxImplicitFunction.h"
#include "MantidGeometry/MDGeometry/MDImplicitFunction.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/CPUTimer.h"
#include "MantidKernel/Memory.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/ProgressText.h"
#include "MantidKernel/ThreadPool.h"
#include "MantidKernel/ThreadScheduler.h"
#include "MantidKernel/Timer.h"
#include "MantidKernel/Utils.h"
#include "MantidAPI/BoxController.h"
#include "MantidDataObjects/CoordTransformDistance.h"
#include "MantidDataObjects/MDBox.h"
#include "MantidDataObjects/MDLeanEvent.h"
#include "MantidDataObjects/MDGridBox.h"
#include <nexus/NeXusFile.hpp>
#include "MantidTestHelpers/MDEventsTestHelper.h"
#include "MDBoxTest.h"
#include <boost/random/linear_congruential.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <cmath>
#include <cxxtest/TestSuite.h>
#include <map>
#include <memory>
#include <Poco/File.h>
#include <vector>
#include <gmock/gmock.h>

using namespace Mantid;
using namespace Mantid::Kernel;
using namespace Mantid::DataObjects;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using namespace testing;

class MDGridBoxTest : public CxxTest::TestSuite {
private:
  /// Mock type to help determine if masking is being determined correctly
  class MockMDBox : public MDBox<MDLeanEvent<1>, 1> {
    API::BoxController *const pBC;

  public:
    MockMDBox()
        : MDBox<MDLeanEvent<1>, 1>(new API::BoxController(1)),
          pBC(MDBox<MDLeanEvent<1>, 1>::getBoxController()) {}
    MOCK_CONST_METHOD0(getIsMasked, bool());
    MOCK_METHOD0(mask, void());
    MOCK_METHOD0(unmask, void());
    ~MockMDBox() { delete pBC; }
  };

  // the sp to a box controller used as general reference to all tested
  // classes/operations including MockMDBox
  BoxController_sptr gbc;

public:
  // This pair of boilerplate methods prevent the suite being created statically
  // This means the constructor isn't called when running other tests
  static MDGridBoxTest *createSuite() { return new MDGridBoxTest(); }
  static void destroySuite(MDGridBoxTest *suite) { delete suite; }

  bool DODEBUG;
  MDGridBoxTest() {
    gbc = BoxController_sptr(new BoxController(1));
    DODEBUG = false;
  }

  //-------------------------------------------------------------------------------------
  void test_MDBoxConstructor() {

    MDBox<MDLeanEvent<1>, 1> *b = MDEventsTestHelper::makeMDBox1(10);
    TS_ASSERT_EQUALS(b->getNumDims(), 1);
    TS_ASSERT_EQUALS(b->getNPoints(), 0);
    TS_ASSERT_DELTA(b->getExtents(0).getMin(), 0.0, 1e-5);
    TS_ASSERT_DELTA(b->getExtents(0).getMax(), 10.0, 1e-5);
    TS_ASSERT_DELTA(b->getVolume(), 10.0, 1e-5);
    // Start at ID 0.
    TS_ASSERT_EQUALS(b->getID(), 0);

    BoxController *const bcc = b->getBoxController();
    delete b;
    if (DODEBUG) {
      std::cout << sizeof(MDLeanEvent<3>) << " bytes per MDLeanEvent(3)"
                << std::endl;
      std::cout << sizeof(MDLeanEvent<4>) << " bytes per MDLeanEvent(4)"
                << std::endl;
      std::cout << sizeof(Mantid::Kernel::Mutex) << " bytes per Mutex"
                << std::endl;
      std::cout << sizeof(MDDimensionExtents<coord_t>)
                << " bytes per MDDimensionExtents" << std::endl;
      std::cout << sizeof(MDBox<MDLeanEvent<3>, 3>) << " bytes per MDBox(3)"
                << std::endl;
      std::cout << sizeof(MDBox<MDLeanEvent<4>, 4>) << " bytes per MDBox(4)"
                << std::endl;
      std::cout << sizeof(MDGridBox<MDLeanEvent<3>, 3>)
                << " bytes per MDGridBox(3)" << std::endl;
      std::cout << sizeof(MDGridBox<MDLeanEvent<4>, 4>)
                << " bytes per MDGridBox(4)" << std::endl;

      MemoryStats mem;
      size_t start = mem.availMem();
      std::cout << start << " KB before" << std::endl;
      CPUTimer tim;
      for (size_t i = 0; i < 1000000; i++) {
        MDBox<MDLeanEvent<3>, 3> *box = new MDBox<MDLeanEvent<3>, 3>(bcc);
        (void)box;
      }
      std::cout << tim << " to allocate a million boxes" << std::endl;
      mem.update();
      size_t stop = mem.availMem();
      std::cout << stop << " KB after " << std::endl;
      std::cout << start - stop << " KB change " << std::endl;
      std::cout << (start - stop) * 1024 / sizeof(MDBox<MDLeanEvent<3>, 3>)
                << " times the sizeof MDBox3" << std::endl;
      delete bcc;
    } else
      delete bcc;
  }

  void check_MDGridBox(MDGridBox<MDLeanEvent<1>, 1> *g) {
    // The grid box stole the ID of the box it replaces.
    TS_ASSERT_EQUALS(g->getID(), 0);

    // Look overall; it has 10 points
    TS_ASSERT_EQUALS(g->getNumDims(), 1);
    TS_ASSERT_EQUALS(g->getNPoints(), 10);
    // Its depth level should be 0 (same as parent)
    TS_ASSERT_EQUALS(g->getDepth(), 0);
    // It was split into 10 MDBoxes.
    TS_ASSERT_EQUALS(g->getNumMDBoxes(), 10);
    // Same result for non-recursive children
    TS_ASSERT_EQUALS(g->getNumChildren(), 10);
    // The volume was set correctly
    TS_ASSERT_DELTA(g->getVolume(), 10.0, 1e-5);

    // It has a BoxController
    TS_ASSERT(g->getBoxController());

    // Check the boxes
    std::vector<MDBoxBase<MDLeanEvent<1>, 1> *> &boxes = g->getBoxes();
    TS_ASSERT_EQUALS(boxes.size(), 10);
    for (size_t i = 0; i < boxes.size(); i++) {
      // The get child method is equivalent
      TS_ASSERT_EQUALS(boxes[i], g->getChild(i));
      MDBox<MDLeanEvent<1>, 1> *box =
          dynamic_cast<MDBox<MDLeanEvent<1>, 1> *>(boxes[i]);

      // Sequential ID, starting at 1 since 0 was used by the parent.
      TS_ASSERT_EQUALS(box->getID(), i + 1);
      // At the right place?
      TS_ASSERT_DELTA(box->getExtents(0).getMin(), double(i) * 1.0, 1e-6);
      TS_ASSERT_DELTA(box->getExtents(0).getMax(), double(i + 1) * 1.0, 1e-6);
      // Look at the single event in there
      TS_ASSERT_EQUALS(box->getNPoints(), 1);
      MDLeanEvent<1> ev = box->getEvents()[0];
      TS_ASSERT_DELTA(ev.getCenter(0), double(i) * 1.0 + 0.5, 1e-5);
      // Its depth level should be 1 (deeper than parent)
      TS_ASSERT_EQUALS(box->getDepth(), 1);
      // The volume was set correctly
      TS_ASSERT_DELTA(box->getVolume(), 1.0, 1e-5);
      // The parent of the MDBox is the grid box
      TS_ASSERT_EQUALS(box->getParent(), g);
    }
  }

  //-------------------------------------------------------------------------------------
  void test_MDGridBox_constructor_from_MDBox() {
    MDBox<MDLeanEvent<1>, 1> *b = MDEventsTestHelper::makeMDBox1();
    TS_ASSERT(b->getBoxController());
    // Start at ID 0.
    TS_ASSERT_EQUALS(b->getID(), 0);
    // Give it 10 events
    const std::vector<MDLeanEvent<1>> events =
        MDEventsTestHelper::makeMDEvents1(10);
    b->addEvents(events);
    TS_ASSERT_EQUALS(b->getNPoints(), 10);
    TS_ASSERT_DELTA(b->getVolume(), 10.0, 1e-5);

    // Build the grid box out of it
    MDGridBox<MDLeanEvent<1>, 1> *g = new MDGridBox<MDLeanEvent<1>, 1>(b);

    // Perform a detailed check
    check_MDGridBox(g);

    // Now we add 10 more events
    // auto events = MDEventsTestHelper::makeMDEvents1(10);
    TSM_ASSERT_EQUALS("No bad events ", 0, g->addEvents(events));

    // And now there should be 2 events per box
    std::vector<MDBoxBase<MDLeanEvent<1>, 1> *> boxes = g->getBoxes();
    for (size_t i = 0; i < 10; i++) {
      MDBox<MDLeanEvent<1>, 1> *box =
          dynamic_cast<MDBox<MDLeanEvent<1>, 1> *>(boxes[i]);
      TS_ASSERT_EQUALS(box->getNPoints(), 2);
    }

    std::vector<signal_t> sigErr(20);
    std::vector<coord_t> coord(10);
    std::vector<uint16_t> runIndex;
    std::vector<uint32_t> detID;

    for (size_t i = 0; i < 10; i++) {
      sigErr[2 * i] = events[i].getSignal();
      sigErr[2 * i + 1] = events[i].getErrorSquared();
      coord[i] = events[i].getCenter(0);
    }

    g->buildAndAddEvents(sigErr, coord, runIndex, detID);

    for (size_t i = 0; i < 10; i++) {
      MDBox<MDLeanEvent<1>, 1> *box =
          dynamic_cast<MDBox<MDLeanEvent<1>, 1> *>(boxes[i]);
      TS_ASSERT_EQUALS(box->getNPoints(), 3);
    }

    BoxController *const bcc = b->getBoxController();
    delete b;
    delete bcc;
    delete g;
  }

  //-------------------------------------------------------------------------------------
  void test_MDGridBox_copy_constructor() {
    MDBox<MDLeanEvent<1>, 1> *b = MDEventsTestHelper::makeMDBox1(10);
    TS_ASSERT_EQUALS(b->getID(), 0);
    const std::vector<MDLeanEvent<1>> events =
        MDEventsTestHelper::makeMDEvents1(10);
    b->addEvents(events);
    TS_ASSERT_EQUALS(b->getNPoints(), 10);
    TS_ASSERT_DELTA(b->getVolume(), 10.0, 1e-5);

    // Build the grid box out of it
    MDGridBox<MDLeanEvent<1>, 1> *g1 = new MDGridBox<MDLeanEvent<1>, 1>(b);
    MDGridBox<MDLeanEvent<1>, 1> *g2 =
        new MDGridBox<MDLeanEvent<1>, 1>(*g1, g1->getBoxController());

    // Perform a detailed check
    check_MDGridBox(g2);

    BoxController *const bcc = b->getBoxController();
    delete bcc;
    delete g2;
    delete g1;
    delete b;
  }

  void test_setBoxController() {
    MDGridBox<MDLeanEvent<1>, 1> *box =
        MDEventsTestHelper::makeMDGridBox<1>(10, 10, 0.0, 10.0);
    BoxController *originalBoxController = box->getBoxController();
    BoxController *const newBoxController = originalBoxController->clone();

    TS_ASSERT_DIFFERS(originalBoxController, newBoxController);

    MDGridBox<MDLeanEvent<1>, 1> *box1 =
        new MDGridBox<MDLeanEvent<1>, 1>(*box, newBoxController);

    auto boxes = box1->getBoxes();
    for (size_t i = 0; i < boxes.size(); ++i) {
      TSM_ASSERT_EQUALS(
          "All child boxes should have the same box controller as the parent.",
          newBoxController, boxes[i]->getBoxController());
    }
    delete newBoxController;
    delete box1;
    delete originalBoxController;
    delete box;
  }

  //-----------------------------------------------------------------------------------------
  /** Manually setting children of a grid box (for NXS file loading) */
  void test_setChildren() {
    // Build the grid box
    MDGridBox<MDLeanEvent<1>, 1> *g =
        MDEventsTestHelper::makeMDGridBox<1>(10, 10, 0.0, 10.0);
    // Clear the initial children
    for (size_t i = 0; i < g->getNumChildren(); i++)
      delete g->getChild(i);

    BoxController *const bcc = g->getBoxController();
    std::vector<API::IMDNode *> boxes;
    for (size_t i = 0; i < 15; i++)
      boxes.push_back(MDEventsTestHelper::makeMDBox1(10, bcc));
    TS_ASSERT_THROWS_NOTHING(g->setChildren(boxes, 2, 12));

    TS_ASSERT_EQUALS(g->getNumChildren(), 10);
    for (size_t i = 2; i < 12; i++) {
      TS_ASSERT_EQUALS(g->getChild(i - 2), boxes[i]);
      // Parent was set correctly in child
      TS_ASSERT_EQUALS(g->getChild(i - 2)->getParent(), g);
    }
    // MDGridBox will delete the children that it pulled in but the rest need to
    // be
    // taken care of manually
    size_t indices[5] = {0, 1, 12, 13, 14};
    for (size_t i = 0; i < 5; ++i)
      delete boxes[indices[i]];
    delete g;
    delete bcc;
  }

  void test_getChildIndexFromID() {
    // Build the grid box
    MDGridBox<MDLeanEvent<1>, 1> *g =
        MDEventsTestHelper::makeMDGridBox<1>(10, 10, 0.0, 10.0);
    TS_ASSERT_EQUALS(g->getChildIndexFromID(g->getChild(0)->getID()), 0);
    TS_ASSERT_EQUALS(g->getChildIndexFromID(g->getChild(5)->getID()), 5);
    TS_ASSERT_EQUALS(g->getChildIndexFromID(0), UNDEF_SIZET);
    TS_ASSERT_EQUALS(g->getChildIndexFromID(11), UNDEF_SIZET);
    BoxController *const bcc = g->getBoxController();
    delete g;
    delete bcc;
  }

  //-------------------------------------------------------------------------------------
  /** Build a 3D MDGridBox and check that the boxes created within are where you
   * expect */
  void test_MDGridBox3() {
    MDBox<MDLeanEvent<3>, 3> *b = MDEventsTestHelper::makeMDBox3();
    // Build the grid box out of it
    MDGridBox<MDLeanEvent<3>, 3> *g = new MDGridBox<MDLeanEvent<3>, 3>(b);
    TS_ASSERT_EQUALS(g->getNumDims(), 3);

    // Check the boxes
    std::vector<MDBoxBase<MDLeanEvent<3>, 3> *> boxes = g->getBoxes();
    TS_ASSERT_EQUALS(boxes.size(), 10 * 5 * 2);
    for (size_t i = 0; i < boxes.size(); i++) {
      MDBox<MDLeanEvent<3>, 3> *box =
          dynamic_cast<MDBox<MDLeanEvent<3>, 3> *>(boxes[i]);
      TS_ASSERT(box);
    }
    MDBox<MDLeanEvent<3>, 3> *box;
    box = dynamic_cast<MDBox<MDLeanEvent<3>, 3> *>(boxes[1]);
    MDEventsTestHelper::extents_match(box, 0, 1.0, 2.0);
    MDEventsTestHelper::extents_match(box, 1, 0.0, 2.0);
    MDEventsTestHelper::extents_match(box, 2, 0.0, 5.0);
    box = dynamic_cast<MDBox<MDLeanEvent<3>, 3> *>(boxes[10]);
    MDEventsTestHelper::extents_match(box, 0, 0.0, 1.0);
    MDEventsTestHelper::extents_match(box, 1, 2.0, 4.0);
    MDEventsTestHelper::extents_match(box, 2, 0.0, 5.0);
    box = dynamic_cast<MDBox<MDLeanEvent<3>, 3> *>(boxes[53]);
    MDEventsTestHelper::extents_match(box, 0, 3.0, 4.0);
    MDEventsTestHelper::extents_match(box, 1, 0.0, 2.0);
    MDEventsTestHelper::extents_match(box, 2, 5.0, 10.0);

    BoxController *const bcc = b->getBoxController();
    delete b;
    delete bcc;
    delete g;
  }

  //-------------------------------------------------------------------------------------
  /** Start with a grid box, split some of its contents into sub-gridded boxes.
   */
  void test_splitContents() {
    MDGridBox<MDLeanEvent<2>, 2> *superbox =
        MDEventsTestHelper::makeMDGridBox<2>();
    MDGridBox<MDLeanEvent<2>, 2> *gb;
    MDBox<MDLeanEvent<2>, 2> *b;

    std::vector<MDBoxBase<MDLeanEvent<2>, 2> *> boxes;

    // Start with 100 boxes
    TS_ASSERT_EQUALS(superbox->getNumMDBoxes(), 100);
    // And ID 0
    TS_ASSERT_EQUALS(superbox->getID(), 0);

    // The box is a MDBox at first
    boxes = superbox->getBoxes();
    b = dynamic_cast<MDBox<MDLeanEvent<2>, 2> *>(boxes[0]);
    TS_ASSERT(b);
    TS_ASSERT_DELTA(b->getVolume(), 1.0, 1e-5);

    // It is the first child, so ID is 1
    TS_ASSERT_EQUALS(b->getID(), 1);
    // There were 101 assigned IDs
    TS_ASSERT_EQUALS(b->getBoxController()->getMaxId(), 100 + 1);

    TS_ASSERT_THROWS_NOTHING(superbox->splitContents(0));

    // Now, it has turned into a GridBox
    boxes = superbox->getBoxes();
    gb = dynamic_cast<MDGridBox<MDLeanEvent<2>, 2> *>(boxes[0]);
    TS_ASSERT(gb);
    TS_ASSERT_DELTA(gb->getVolume(), 1.0, 1e-5);

    // ID of first child remains unchanged at 1
    TS_ASSERT_EQUALS(gb->getID(), 1);
    // There were 101 assigned IDs
    TS_ASSERT_EQUALS(gb->getBoxController()->getMaxId(), 200 + 1);
    // The first child of the sub-divided box got 101 as its id
    TS_ASSERT_EQUALS(gb->getBoxes()[0]->getID(), 101);

    // There are now 199 MDBoxes; the 99 at level 1, and 100 at level 2
    TS_ASSERT_EQUALS(superbox->getNumMDBoxes(), 199);

    // You can split it again and it does nothing
    TS_ASSERT_THROWS_NOTHING(superbox->splitContents(0));

    // Still a grid box
    boxes = superbox->getBoxes();
    gb = dynamic_cast<MDGridBox<MDLeanEvent<2>, 2> *>(boxes[0]);
    TS_ASSERT(gb);

    BoxController *const bcc = superbox->getBoxController();
    delete superbox;
    delete bcc;
  }

  //-------------------------------------------------------------------------------------
  /** Adding a single event pushes it as deep as the current grid
   * hierarchy allows
   */
  void test_addEvent_with_recursive_gridding() {
    MDGridBox<MDLeanEvent<2>, 2> *gb;
    MDBox<MDLeanEvent<2>, 2> *b;
    std::vector<MDBoxBase<MDLeanEvent<2>, 2> *> boxes;

    // 10x10 box, extents 0-10.0
    MDGridBox<MDLeanEvent<2>, 2> *superbox =
        MDEventsTestHelper::makeMDGridBox<2>();
    // And the 0-th box is further split (
    TS_ASSERT_THROWS_NOTHING(superbox->splitContents(0));

    TS_ASSERT_EQUALS(superbox->getNPoints(), 0);
    { // One event in 0th box of the 0th box.
      double centers[2] = {0.05, 0.05};
      superbox->addEvent(MDLeanEvent<2>(2.0, 2.0, centers));
    }
    { // One event in 1st box of the 0th box.
      double centers[2] = {0.15, 0.05};
      superbox->addEvent(MDLeanEvent<2>(2.0, 2.0, centers));
    }
    { // One event in 99th box.
      double centers[2] = {9.5, 9.5};
      superbox->addEvent(MDLeanEvent<2>(2.0, 2.0, centers));
    }

    // You must refresh the cache after adding individual events.
    superbox->refreshCache(NULL);
    // superbox->refreshCentroid(NULL);

    TS_ASSERT_EQUALS(superbox->getNPoints(), 3);

    std::vector<coord_t> cenroid(2, 0);
    TS_ASSERT_THROWS(superbox->calculateCentroid(&cenroid[0]),
                     std::runtime_error);

    // Check the centroid for these 3 events
    // TS_ASSERT_DELTA( cenroid[0], 3.233, 0.001);
    // TS_ASSERT_DELTA(cenroid[1] , 3.200, 0.001);

    { // One event in 0th box of the 0th box.
      std::vector<coord_t> centers(2, 0.05f);
      superbox->buildAndAddEvent(2., 2., centers, 0, 0);
    }
    { // One event in 1st box of the 0th box.
      std::vector<coord_t> centers(2, 0.05f);
      centers[0] = 0.15f;
      superbox->buildAndAddEvent(2., 2., centers, 0, 0);
    }
    { // One event in 99th box.
      std::vector<coord_t> centers(2, 9.5);
      superbox->buildAndAddEvent(2.0, 2.0, centers, 0, 0);
    }
    TS_ASSERT_EQUALS(superbox->getNPoints(), 3);

    superbox->refreshCache(NULL);
    TS_ASSERT_EQUALS(superbox->getNPoints(), 6);

    // Retrieve the 0th grid box
    boxes = superbox->getBoxes();
    gb = dynamic_cast<MDGridBox<MDLeanEvent<2>, 2> *>(boxes[0]);
    TS_ASSERT(gb);

    // It has three points
    TS_ASSERT_EQUALS(gb->getNPoints(), 4);

    // Retrieve the MDBox at 0th and 1st indexes in THAT gridbox
    boxes = gb->getBoxes();
    b = dynamic_cast<MDBox<MDLeanEvent<2>, 2> *>(boxes[0]);
    TS_ASSERT_EQUALS(b->getNPoints(), 2);
    b = dynamic_cast<MDBox<MDLeanEvent<2>, 2> *>(boxes[1]);
    TS_ASSERT_EQUALS(b->getNPoints(), 2);

    // Get the 99th box at the first level. It is not split
    boxes = superbox->getBoxes();
    b = dynamic_cast<MDBox<MDLeanEvent<2>, 2> *>(boxes[99]);
    TS_ASSERT(b);
    if (!b)
      return;
    // And it has only the one point
    TS_ASSERT_EQUALS(b->getNPoints(), 2);

    BoxController *const bcc = superbox->getBoxController();
    delete superbox;
    delete bcc;
  }

  //-------------------------------------------------------------------------------------
  void test_transformDimensions() {
    MDBox<MDLeanEvent<1>, 1> *b = MDEventsTestHelper::makeMDBox1();

    // Give it 10 events
    const std::vector<MDLeanEvent<1>> events =
        MDEventsTestHelper::makeMDEvents1(10);
    b->addEvents(events);
    MDGridBox<MDLeanEvent<1>, 1> *g = new MDGridBox<MDLeanEvent<1>, 1>(b);
    TSM_ASSERT_EQUALS("MDBoxes start with 1 each.",
                      g->getChild(9)->getNPoints(), 1);

    std::vector<double> scaling(1, 3.0);
    std::vector<double> offset(1, 1.0);
    g->transformDimensions(scaling, offset);

    TS_ASSERT_DELTA(g->getVolume(), 30.0, 1e-5);
    MDLeanEvent<1> ev;
    ev.setCenter(0, 30.9f);
    g->addEvent(ev);
    TSM_ASSERT_EQUALS("New event was added in the right spot.",
                      g->getChild(9)->getNPoints(), 2);

    BoxController *const bcc = b->getBoxController();
    delete b;
    delete bcc;
    delete g;
  }

  //-------------------------------------------------------------------------------------
  /** Recursive getting of a list of MDBoxBase */
  void test_getBoxes() {
    MDGridBox<MDLeanEvent<1>, 1> *parent =
        MDEventsTestHelper::makeRecursiveMDGridBox<1>(3, 3);
    TS_ASSERT(parent);
    std::vector<API::IMDNode *> boxes;

    boxes.clear();
    parent->getBoxes(boxes, 0, false);
    TS_ASSERT_EQUALS(boxes.size(), 1);
    TS_ASSERT_EQUALS(boxes[0], parent);

    boxes.clear();
    parent->getBoxes(boxes, 1, false);
    TS_ASSERT_EQUALS(boxes.size(), 4);
    TS_ASSERT_EQUALS(boxes[0], parent);
    TS_ASSERT_EQUALS(boxes[1]->getDepth(), 1);

    boxes.clear();
    parent->getBoxes(boxes, 2, false);
    TS_ASSERT_EQUALS(boxes.size(), 4 + 9);
    TS_ASSERT_EQUALS(boxes[0], parent);
    TS_ASSERT_EQUALS(boxes[1]->getDepth(), 1);
    TS_ASSERT_EQUALS(boxes[2]->getDepth(), 2);

    boxes.clear();
    parent->getBoxes(boxes, 3, false);
    TS_ASSERT_EQUALS(boxes.size(), 4 + 9 + 27);

    // Leaves only = only report the deepest boxes.
    boxes.clear();
    parent->getBoxes(boxes, 3, true);
    TS_ASSERT_EQUALS(boxes.size(), 27);
    TS_ASSERT_EQUALS(boxes[0]->getDepth(), 3);

    // Leaves only, with limited depth = report the max depth if that is the
    // effective 'leaf'
    boxes.clear();
    parent->getBoxes(boxes, 2, true);
    TS_ASSERT_EQUALS(boxes.size(), 9);
    TS_ASSERT_EQUALS(boxes[0]->getDepth(), 2);

    BoxController *const bcc = parent->getBoxController();
    delete parent;
    delete bcc;
  }

  //-------------------------------------------------------------------------------------
  /** Recursive getting of a list of MDBoxBase, with an implicit function
   * limiting it */
  void test_getBoxes_ImplicitFunction() {

    MDGridBox<MDLeanEvent<1>, 1> *parent =
        MDEventsTestHelper::makeRecursiveMDGridBox<1>(4, 3);
    TS_ASSERT(parent);
    std::vector<API::IMDNode *> boxes;

    // Function of everything x > 1.51
    MDImplicitFunction *function = new MDImplicitFunction;
    coord_t normal[1] = {+1};
    coord_t origin[1] = {1.51f};
    function->addPlane(MDPlane(1, normal, origin));

    boxes.clear();
    parent->getBoxes(boxes, 3, false, function);
    TS_ASSERT_EQUALS(boxes.size(), 54);
    // The boxes extents make sense
    for (size_t i = 0; i < boxes.size(); i++) {
      TS_ASSERT(boxes[i]->getExtents(0).getMax() >= 1.51);
    }

    // --- Now leaf-only ---
    boxes.clear();
    parent->getBoxes(boxes, 3, true, function);
    TS_ASSERT_EQUALS(boxes.size(), 40);
    // The boxes extents make sense
    for (size_t i = 0; i < boxes.size(); i++) {
      TS_ASSERT(boxes[i]->getExtents(0).getMax() >= 1.51);
    }

    // Limit by another plane
    coord_t normal2[1] = {-1};
    coord_t origin2[1] = {2.99f};
    function->addPlane(MDPlane(1, normal2, origin2));
    boxes.clear();
    parent->getBoxes(boxes, 3, false, function);
    TS_ASSERT_EQUALS(boxes.size(), 33);
    for (size_t i = 0; i < boxes.size(); i++) {
      TS_ASSERT(boxes[i]->getExtents(0).getMax() >= 1.51);
      TS_ASSERT(boxes[i]->getExtents(0).getMin() <= 2.99);
    }

    // Same, leaf only
    boxes.clear();
    parent->getBoxes(boxes, 3, true, function);
    TS_ASSERT_EQUALS(boxes.size(), 24);
    for (size_t i = 0; i < boxes.size(); i++) {
      TS_ASSERT(boxes[i]->getExtents(0).getMax() >= 1.51);
      TS_ASSERT(boxes[i]->getExtents(0).getMin() <= 2.99);
    }

    // ----- Infinitely thin plane for an implicit function ------------
    delete function;
    function = new MDImplicitFunction;
    coord_t normal3[1] = {-1};
    coord_t origin3[1] = {1.51f};
    function->addPlane(MDPlane(1, normal, origin));
    function->addPlane(MDPlane(1, normal3, origin3));

    boxes.clear();
    parent->getBoxes(boxes, 3, true, function);
    TSM_ASSERT_EQUALS("Only one box is found by an infinitely thin plane",
                      boxes.size(), 1);

    // clean up  behind
    BoxController *const bcc = parent->getBoxController();
    delete parent;
    delete bcc;
    delete function;
  }

  //-------------------------------------------------------------------------------------
  /** Recursive getting of a list of MDBoxBase, with an implicit function
   * limiting it */
  void test_getBoxes_ImplicitFunction_2D() {

    MDGridBox<MDLeanEvent<2>, 2> *parent =
        MDEventsTestHelper::makeRecursiveMDGridBox<2>(4, 1);
    TS_ASSERT(parent);
    std::vector<API::IMDNode *> boxes;

    // Function of x,y between 2 and 3
    std::vector<coord_t> min(2, 1.99f);
    std::vector<coord_t> max(2, 3.01f);
    MDImplicitFunction *function = new MDBoxImplicitFunction(min, max);

    boxes.clear();
    parent->getBoxes(boxes, 3, false, function);
    TS_ASSERT_EQUALS(boxes.size(), 46);
    // The boxes extents make sense
    for (size_t i = 0; i < boxes.size(); i++) {
      TS_ASSERT(boxes[i]->getExtents(0).getMax() >= 2.00);
      TS_ASSERT(boxes[i]->getExtents(0).getMin() <= 3.00);
      TS_ASSERT(boxes[i]->getExtents(1).getMax() >= 2.00);
      TS_ASSERT(boxes[i]->getExtents(1).getMin() <= 3.00);
    }

    // -- Leaf only ---
    boxes.clear();
    parent->getBoxes(boxes, 3, true, function);
    TS_ASSERT_EQUALS(
        boxes.size(),
        16 + 4 * 4 +
            4); // 16 in the center one + 4x4 at the 4 edges + 4 at the corners
    // The boxes extents make sense
    for (size_t i = 0; i < boxes.size(); i++) {
      TS_ASSERT(boxes[i]->getExtents(0).getMax() >= 2.00);
      TS_ASSERT(boxes[i]->getExtents(0).getMin() <= 3.00);
      TS_ASSERT(boxes[i]->getExtents(1).getMax() >= 2.00);
      TS_ASSERT(boxes[i]->getExtents(1).getMin() <= 3.00);
    }

    // clean up  behind
    BoxController *const bcc = parent->getBoxController();
    delete parent;
    delete bcc;
    delete function;
  }

  //-------------------------------------------------------------------------------------
  /** Recursive getting of a list of MDBoxBase, with an implicit function
   * limiting it.
   * This implicit function is a box of size 0. */
  void test_getBoxes_ZeroSizeImplicitFunction_2D() {
    MDGridBox<MDLeanEvent<2>, 2> *parent =
        MDEventsTestHelper::makeRecursiveMDGridBox<2>(4, 1);
    TS_ASSERT(parent);
    std::vector<API::IMDNode *> boxes;

    // Function of x,y with 0 width and height
    std::vector<coord_t> min(2, 1.99f);
    std::vector<coord_t> max(2, 1.99f);
    MDImplicitFunction *function = new MDBoxImplicitFunction(min, max);

    boxes.clear();
    parent->getBoxes(boxes, 3, false, function);
    // Returns 3 boxes: the big one with everything, the size 1 under, and the
    // size 0.25 under that
    TS_ASSERT_EQUALS(boxes.size(), 3);

    // Leaf only: returns only one box
    boxes.clear();
    parent->getBoxes(boxes, 3, true, function);
    TS_ASSERT_EQUALS(boxes.size(), 1);
    TS_ASSERT_DELTA(boxes[0]->getExtents(0).getMin(), 1.75, 1e-4);
    TS_ASSERT_DELTA(boxes[0]->getExtents(0).getMax(), 2.00, 1e-4);

    // clean up  behind
    BoxController *const bcc = parent->getBoxController();
    delete parent;
    delete bcc;
    delete function;
  }

  //-------------------------------------------------------------------------------------
  /** Recursive getting of a list of MDBoxBase, with an implicit function
   * limiting it.
   * This implicit function is a box of size 0. */
  void test_getBoxes_ZeroSizeImplicitFunction_4D() {
    MDGridBox<MDLeanEvent<4>, 4> *parent =
        MDEventsTestHelper::makeRecursiveMDGridBox<4>(4, 1);
    TS_ASSERT(parent);
    std::vector<API::IMDNode *> boxes;

    // Function of x,y with 0 width and height
    std::vector<coord_t> min(4, 1.99f);
    std::vector<coord_t> max(4, 1.99f);
    MDImplicitFunction *function = new MDBoxImplicitFunction(min, max);

    boxes.clear();
    parent->getBoxes(boxes, 3, false, function);
    // Returns 3 boxes: the big one with everything, the size 1 under, and the
    // size 0.25 under that
    TS_ASSERT_EQUALS(boxes.size(), 3);

    // Leaf only: returns only one box
    boxes.clear();
    parent->getBoxes(boxes, 3, true, function);
    TS_ASSERT_EQUALS(boxes.size(), 1);
    TS_ASSERT_DELTA(boxes[0]->getExtents(0).getMin(), 1.75, 1e-4);
    TS_ASSERT_DELTA(boxes[0]->getExtents(0).getMax(), 2.00, 1e-4);

    // clean up  behind
    BoxController *const bcc = parent->getBoxController();
    delete parent;
    delete bcc;
    delete function;
  }

  //-------------------------------------------------------------------------------------
  /** Gauge how fast addEvent is with several levels of gridding
   * NOTE: DISABLED because it is slow.
   * */
  void xtest_addEvent_with_recursive_gridding_Performance() {
    // Make a 2D box split into 4, 4 levels deep. = 4^4^2 boxes at the bottom =
    // 256^2 boxes.
    size_t numSplit = 4;
    for (size_t recurseLevels = 1; recurseLevels < 5; recurseLevels++) {
      double boxes_per_side = pow(double(numSplit), double(recurseLevels));
      double spacing = double(numSplit) / boxes_per_side;
      // How many times to add the same event
      size_t num_to_repeat = size_t(1e7 / (boxes_per_side * boxes_per_side));

      MDGridBox<MDLeanEvent<2>, 2> *box =
          MDEventsTestHelper::makeRecursiveMDGridBox<2>(numSplit,
                                                        recurseLevels);

      for (double x = 0; x < numSplit; x += spacing)
        for (double y = 0; y < numSplit; y += spacing) {
          for (size_t i = 0; i < num_to_repeat; i++) {
            coord_t centers[2] = {static_cast<coord_t>(x),
                                  static_cast<coord_t>(y)};
            box->addEvent(MDLeanEvent<2>(2.0, 2.0, centers));
          }
        }
      // You must refresh the cache after adding individual events.
      box->refreshCache(NULL);
    }
  }

  //-------------------------------------------------------------------------------------
  /** Fill a 10x10 gridbox with events
   *
   * Tests that bad events are thrown out when using addEvents.
   * */
  void test_addEvents_2D() {
    MDGridBox<MDLeanEvent<2>, 2> *b = MDEventsTestHelper::makeMDGridBox<2>();
    std::vector<MDLeanEvent<2>> events;

    // Make an event in the middle of each box
    for (double x = 0.5; x < 10; x += 1.0)
      for (double y = 0.5; y < 10; y += 1.0) {
        coord_t centers[2] = {static_cast<coord_t>(x), static_cast<coord_t>(y)};
        events.push_back(MDLeanEvent<2>(2.0, 2.0, centers));
      }

    size_t numbad = 0;
    TS_ASSERT_THROWS_NOTHING(numbad = b->addEvents(events););
    // Get the right totals again
    b->refreshCache(NULL);
    TS_ASSERT_EQUALS(numbad, 0);
    TS_ASSERT_EQUALS(b->getNPoints(), 100);
    TS_ASSERT_EQUALS(b->getSignal(), 100 * 2.0);
    TS_ASSERT_EQUALS(b->getErrorSquared(), 100 * 2.0);
    TS_ASSERT_DELTA(b->getSignalNormalized(), 100 * 2.0 / 100.0, 1e-5);
    TS_ASSERT_DELTA(b->getErrorSquaredNormalized(), 100 * 2.0 / 100.0, 1e-5);

    // Get all the boxes contained
    std::vector<MDBoxBase<MDLeanEvent<2>, 2> *> boxes = b->getBoxes();
    TS_ASSERT_EQUALS(boxes.size(), 100);
    for (size_t i = 0; i < boxes.size(); i++) {
      TS_ASSERT_EQUALS(boxes[i]->getNPoints(), 1);
      TS_ASSERT_EQUALS(boxes[i]->getSignal(), 2.0);
      TS_ASSERT_EQUALS(boxes[i]->getErrorSquared(), 2.0);
      TS_ASSERT_EQUALS(boxes[i]->getSignalNormalized(), 2.0);
      TS_ASSERT_EQUALS(boxes[i]->getErrorSquaredNormalized(), 2.0);
    }

    // Now try to add bad events (outside bounds)
    events.clear();
    for (double x = -5.0; x < 20; x += 20.0)
      for (double y = -5.0; y < 20; y += 20.0) {
        double centers[2] = {x, y};
        events.push_back(MDLeanEvent<2>(2.0, 2.0, centers));
      }
    // Get the right totals again
    b->refreshCache(NULL);
    // All 4 points get rejected
    TS_ASSERT_THROWS_NOTHING(numbad = b->addEvents(events););
    TS_ASSERT_EQUALS(numbad, 4);
    // Number of points and signal is unchanged.
    TS_ASSERT_EQUALS(b->getNPoints(), 100);
    TS_ASSERT_EQUALS(b->getSignal(), 100 * 2.0);
    TS_ASSERT_EQUALS(b->getErrorSquared(), 100 * 2.0);

    // clean up  behind
    BoxController *const bcc = b->getBoxController();
    delete b;
    delete bcc;
  }

  ////-------------------------------------------------------------------------------------
  ///** Tests add_events with limits into the vectorthat bad events are thrown
  /// out when using addEvents.
  // * */
  // void xest_addEvents_start_stop()
  //{
  //  MDGridBox<MDLeanEvent<2>,2> * b = MDEventsTestHelper::makeMDGridBox<2>();
  //  std::vector< MDLeanEvent<2> > events;

  //  // Make an event in the middle of each box
  //  for (double x=0.5; x < 10; x += 1.0)
  //    for (double y=0.5; y < 10; y += 1.0)
  //    {
  //      double centers[2] = {x,y};
  //      events.push_back( MDLeanEvent<2>(2.0, 2.0, centers) );
  //    }

  //  size_t numbad = 0;
  //  TS_ASSERT_THROWS_NOTHING( numbad = b->addEventsPart( events, 50, 60 ); );
  //  // Get the right totals again
  //  b->refreshCache(NULL);
  //  TS_ASSERT_EQUALS( numbad, 0);
  //  TS_ASSERT_EQUALS( b->getNPoints(), 10);
  //  TS_ASSERT_EQUALS( b->getSignal(), 10*2.0);
  //  TS_ASSERT_EQUALS( b->getErrorSquared(), 10*2.0);
  //}

  //-------------------------------------------------------------------------------------
  /** Test that adding events (as vectors) in parallel does not cause
   * segfaults or incorrect totals.
   * */
  void do_test_addEvents_inParallel(ThreadScheduler *ts) {
    MDGridBox<MDLeanEvent<2>, 2> *b = MDEventsTestHelper::makeMDGridBox<2>();
    int num_repeat = 1000;

    PARALLEL_FOR_NO_WSP_CHECK()
    for (int i = 0; i < num_repeat; i++) {
      std::vector<MDLeanEvent<2>> events;
      // Make an event in the middle of each box
      for (double x = 0.5; x < 10; x += 1.0)
        for (double y = 0.5; y < 10; y += 1.0) {
          double centers[2] = {x, y};
          events.push_back(MDLeanEvent<2>(2.0, 2.0, centers));
        }
      TS_ASSERT_THROWS_NOTHING(b->addEvents(events););
    }

    // Get the right totals again by refreshing
    b->refreshCache(ts);
    TS_ASSERT_EQUALS(b->getNPoints(), 100 * num_repeat);
    TS_ASSERT_EQUALS(b->getSignal(), 100 * num_repeat * 2.0);
    TS_ASSERT_EQUALS(b->getErrorSquared(), 100 * num_repeat * 2.0);

    // clean up  behind
    BoxController *const bcc = b->getBoxController();
    delete b;
    delete bcc;
  }

  void test_addEvents_inParallel() { do_test_addEvents_inParallel(NULL); }

  /** Disabled because parallel RefreshCache is not implemented. Might not be
   * ever? */
  void xtest_addEvents_inParallel_then_refreshCache_inParallel() {
    ThreadScheduler *ts = new ThreadSchedulerFIFO();
    do_test_addEvents_inParallel(ts);
    ThreadPool tp(ts);
    tp.joinAll();
  }

  //-------------------------------------------------------------------------------------
  /** Get a sub-box at a given coord */
  void test_getBoxAtCoord() {
    MDGridBox<MDLeanEvent<2>, 2> *b = MDEventsTestHelper::makeMDGridBox<2>();
    coord_t coords[2] = {1.5, 1.5};
    const MDBoxBase<MDLeanEvent<2>, 2> *c =
        dynamic_cast<const MDBoxBase<MDLeanEvent<2>, 2> *>(
            b->getBoxAtCoord(coords));
    TS_ASSERT_EQUALS(c, b->getChild(11));
    delete b->getBoxController();
    delete b;
  }

  //-------------------------------------------------------------------------------------
  /** Test the routine that auto-splits MDBoxes into MDGridBoxes recursively.
   * It tests the max_depth of splitting too, because there are numerous
   * repeated events at exactly the same position = impossible to separate
   * further.
   * */
  void test_splitAllIfNeeded() {
    typedef MDGridBox<MDLeanEvent<2>, 2> gbox_t;
    typedef MDBox<MDLeanEvent<2>, 2> box_t;
    typedef MDBoxBase<MDLeanEvent<2>, 2> ibox_t;

    gbox_t *b0 = MDEventsTestHelper::makeMDGridBox<2>();
    b0->getBoxController()->setSplitThreshold(100);
    b0->getBoxController()->setMaxDepth(4);

    // Make a 1000 events at exactly the same point
    size_t num_repeat = 1000;
    std::vector<MDLeanEvent<2>> events;
    for (size_t i = 0; i < num_repeat; i++) {
      // Make an event in the middle of each box
      double centers[2] = {1e-10, 1e-10};
      events.push_back(MDLeanEvent<2>(2.0, 2.0, centers));
    }
    TS_ASSERT_THROWS_NOTHING(b0->addEvents(events););

    // Split into sub-grid boxes
    TS_ASSERT_THROWS_NOTHING(b0->splitAllIfNeeded(NULL);)

    // Dig recursively into the gridded box hierarchies
    std::vector<ibox_t *> boxes;
    size_t expected_depth = 0;
    gbox_t *b = b0;
    while (b) {
      expected_depth++;
      boxes = b->getBoxes();

      // Get the 0th box
      b = dynamic_cast<gbox_t *>(boxes[0]);

      // The 0-th box is a MDGridBox (it was split)
      // (though it is normal for b to be a MDBox when you reach the max depth)
      if (expected_depth < 4) {
        TS_ASSERT(b)
      }

      // The 0-th box has all the points
      TS_ASSERT_EQUALS(boxes[0]->getNPoints(), num_repeat);
      // The 0-th box is at the expected_depth
      TS_ASSERT_EQUALS(boxes[0]->getDepth(), expected_depth);