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 "MantidMDEvents/CoordTransformDistance.h"
#include "MantidMDEvents/MDBox.h"
#include "MantidMDEvents/MDLeanEvent.h"
#include "MantidMDEvents/MDGridBox.h"
#include "MantidNexusCPP/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 <cxxtest/TestSuite.h>
#include <iomanip>
#include <map>
#include <memory>
#include <Poco/File.h>
#include <vector>

using namespace Mantid;
using namespace Mantid::Kernel;
using namespace Mantid::MDEvents;
using namespace Mantid::API;
using namespace Mantid::Geometry;

class MDGridBoxTest :    public CxxTest::TestSuite
{
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()
  {
    DODEBUG = false;
  }

  //-------------------------------------------------------------------------------------
  void test_MDBoxConstructor()
  {
    MDBox<MDLeanEvent<1>,1> * b = MDEventsTestHelper::makeMDBox1();
    TS_ASSERT_EQUALS( b->getNumDims(), 1);
    TS_ASSERT_EQUALS( b->getNPoints(), 0);
    TS_ASSERT_DELTA( b->getExtents(0).min, 0.0,  1e-5);
    TS_ASSERT_DELTA( b->getExtents(0).max, 10.0, 1e-5);
    TS_ASSERT_DELTA( b->getVolume(), 10.0, 1e-5);
    // Start at ID 0.
    TS_ASSERT_EQUALS( b->getId(), 0);
    delete b;


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


  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<IMDBox<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).min, double(i)*1.0, 1e-6);
      TS_ASSERT_DELTA(box->getExtents(0).max, 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();
    // 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
    g->addEvents( MDEventsTestHelper::makeMDEvents1(10) );

    // And now there should be 2 events per box
    std::vector<IMDBox<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);
    }
  }


  //-------------------------------------------------------------------------------------
  void test_MDGridBox_copy_constructor()
  {
    MDBox<MDLeanEvent<1>,1> * b = MDEventsTestHelper::makeMDBox1();
    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);

    // Perform a detailed check
    check_MDGridBox(g2);
  }



  //-----------------------------------------------------------------------------------------
  /** Test splitting of a MDBox into a MDGridBox when the
   * original box is backed by a file. */
  void test_fileBackEnd_construction()
  {
    // Create a box with a controller for the back-end
    BoxController_sptr bc(new BoxController(3));
    bc->setSplitInto(5);
    // Handle the disk MRU values
    bc->setCacheParameters(sizeof(MDLeanEvent<3>), 10000);
    DiskBuffer & dbuf = bc->getDiskBuffer();
    // Make a box from 0-10 in 3D
    MDBox<MDLeanEvent<3>,3> * c = new MDBox<MDLeanEvent<3>,3>(bc, 0);
    for (size_t d=0; d<3; d++) c->setExtents(d, 0, 10);

    // Create and open the test NXS file
    ::NeXus::File * file = MDBoxTest::do_saveAndOpenNexus(*c, "MDGridBoxTest.nxs");
    TSM_ASSERT_EQUALS( "1000 events (on file)", c->getNPoints(), 1000);

    // At this point the MDBox is set to be on disk
    TSM_ASSERT_EQUALS( "No free blocks to start with", dbuf.getFreeSpaceMap().size(), 0);

    // Construct the grid box by splitting the MDBox
    MDGridBox<MDLeanEvent<3>,3> * gb = new MDGridBox<MDLeanEvent<3>,3>(c);
    TSM_ASSERT_EQUALS( "Grid box also has 1000 points", gb->getNPoints(), 1000);
    TSM_ASSERT_EQUALS( "Grid box has 125 children (5x5x5)", gb->getNumChildren(), 125);
    TSM_ASSERT_EQUALS( "The old spot in the file is now free", dbuf.getFreeSpaceMap().size(), 1);

    // Get a child
    MDBox<MDLeanEvent<3>,3> * b = dynamic_cast<MDBox<MDLeanEvent<3>,3> *>(gb->getChild(22));
    TSM_ASSERT_EQUALS( "Child has 8 events", b->getNPoints(), 8);
    TSM_ASSERT_EQUALS( "Child is NOT on disk", b->getOnDisk(), false);

    file->close();
    MDBoxTest::do_deleteNexusFile("MDGridBoxTest.nxs");
  }


  //-----------------------------------------------------------------------------------------
  /** 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);
    std::vector<IMDBox<MDLeanEvent<1>,1>*> boxes;
    for (size_t i=0; i<15; i++)
      boxes.push_back( MDEventsTestHelper::makeMDBox1() );
    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);
    }
  }

  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), size_t(-1) );
    TS_ASSERT_EQUALS(g->getChildIndexFromID(11), size_t(-1) );
  }



  //-------------------------------------------------------------------------------------
  /** 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<IMDBox<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);
  }


  //-------------------------------------------------------------------------------------
  /** 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<IMDBox<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 );
  }

  //-------------------------------------------------------------------------------------
  /** 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<IMDBox<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.
      coord_t centers[2] = {0.05, 0.05};
      superbox->addEvent( MDLeanEvent<2>(2.0, 2.0, centers) );
    }
    { // One event in 1st box of the 0th box.
      coord_t centers[2] = {0.15, 0.05};
      superbox->addEvent( MDLeanEvent<2>(2.0, 2.0, centers) );
    }
    { // One event in 99th box.
      coord_t 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 );

#ifdef MDBOX_TRACK_CENTROID
    // Check the centroid for these 3 events
    TS_ASSERT_DELTA( superbox->getCentroid(0), 3.233, 0.001);
    TS_ASSERT_DELTA( superbox->getCentroid(1), 3.200, 0.001);
#endif

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

    // It has two points
    TS_ASSERT_EQUALS( gb->getNPoints(), 2 );

    // 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(), 1 );
    b = dynamic_cast<MDBox<MDLeanEvent<2>,2> *>(boxes[1]);
    TS_ASSERT_EQUALS( b->getNPoints(), 1 );

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

  //-------------------------------------------------------------------------------------
  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.9);
    g->addEvent(ev);
    TSM_ASSERT_EQUALS("New event was added in the right spot.", g->getChild(9)->getNPoints(), 2);
  }



  //-------------------------------------------------------------------------------------
  /** Recursive getting of a list of IMDBox */
  void test_getBoxes()
  {
    MDGridBox<MDLeanEvent<1>,1> * parent = MDEventsTestHelper::makeRecursiveMDGridBox<1>(3,3);
    TS_ASSERT(parent);
    std::vector<IMDBox<MDLeanEvent<1>,1> *> 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);
  }


  //-------------------------------------------------------------------------------------
  /** Recursive getting of a list of IMDBox, 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<IMDBox<MDLeanEvent<1>,1> *> boxes;

    // Function of everything x > 1.51
    MDImplicitFunction * function = new MDImplicitFunction;
    coord_t normal[1] = {+1};
    coord_t origin[1] = {1.51};
    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).max >= 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).max >= 1.51);
    }

    // Limit by another plane
    coord_t normal2[1] = {-1};
    coord_t origin2[1] = {2.99};
    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).max >= 1.51);
      TS_ASSERT( boxes[i]->getExtents(0).min <= 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).max >= 1.51);
      TS_ASSERT( boxes[i]->getExtents(0).min <= 2.99);
    }

    // ----- Infinitely thin plane for an implicit function ------------
    function = new MDImplicitFunction;
    coord_t normal3[1] = {-1};
    coord_t origin3[1] = {1.51};
    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);

  }



  //-------------------------------------------------------------------------------------
  /** Recursive getting of a list of IMDBox, 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<IMDBox<MDLeanEvent<2>,2> *> boxes;

    // Function of x,y between 2 and 3
    std::vector<coord_t> min(2, 1.99);
    std::vector<coord_t> max(2, 3.01);
    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++)
    {
      //std::cout << boxes[i]->getExtentsStr() << std::endl;
      TS_ASSERT( boxes[i]->getExtents(0).max >= 2.00);
      TS_ASSERT( boxes[i]->getExtents(0).min <= 3.00);
      TS_ASSERT( boxes[i]->getExtents(1).max >= 2.00);
      TS_ASSERT( boxes[i]->getExtents(1).min <= 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++)
    {
      //std::cout << boxes[i]->getExtentsStr() << std::endl;
      TS_ASSERT( boxes[i]->getExtents(0).max >= 2.00);
      TS_ASSERT( boxes[i]->getExtents(0).min <= 3.00);
      TS_ASSERT( boxes[i]->getExtents(1).max >= 2.00);
      TS_ASSERT( boxes[i]->getExtents(1).min <= 3.00);
    }

  }



  //-------------------------------------------------------------------------------------
  /** Recursive getting of a list of IMDBox, 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<IMDBox<MDLeanEvent<2>,2> *> boxes;

    // Function of x,y with 0 width and height
    std::vector<coord_t> min(2, 1.99);
    std::vector<coord_t> max(2, 1.99);
    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).min, 1.75, 1e-4);
    TS_ASSERT_DELTA( boxes[0]->getExtents(0).max, 2.00, 1e-4);
  }

  //-------------------------------------------------------------------------------------
  /** Recursive getting of a list of IMDBox, 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<IMDBox<MDLeanEvent<4>,4> *> boxes;

    // Function of x,y with 0 width and height
    std::vector<coord_t> min(4, 1.99);
    std::vector<coord_t> max(4, 1.99);
    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).min, 1.75, 1e-4);
    TS_ASSERT_DELTA( boxes[0]->getExtents(0).max, 2.00, 1e-4);
  }



  //-------------------------------------------------------------------------------------
  /** 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++)
    {
      std::cout << " --- Recursion Level " << recurseLevels << " --- " << std::endl;
      Timer tim1;
      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);
      std::cout << tim1.elapsed() << " seconds to generate the " << boxes_per_side << "^2 boxes." << std::endl;

      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] = {x,y};
            box->addEvent( MDLeanEvent<2>(2.0, 2.0, centers) );
          }
        }
      // You must refresh the cache after adding individual events.
      box->refreshCache(NULL);

      double sec = tim1.elapsed();
      std::cout << sec << " seconds to add " << box->getNPoints() << " events. Each box had " << num_to_repeat << " events." << std::endl;
      std::cout << "equals " << 1e6*sec/double(box->getNPoints()) << " seconds per million events." << std::endl;
    }

  }



  //-------------------------------------------------------------------------------------
  /** 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] = {x,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<IMDBox<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)
      {
        coord_t 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);
  }


  //-------------------------------------------------------------------------------------
  /** Tests add_events with limits into the vectorthat bad events are thrown out when using addEvents.
   * */
  void test_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)
      {
        coord_t centers[2] = {x,y};
        events.push_back( MDLeanEvent<2>(2.0, 2.0, centers) );
      }

    size_t numbad = 0;
    TS_ASSERT_THROWS_NOTHING( numbad = b->addEvents( 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)
        {
          coord_t 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);
  }


  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 IMDBox<MDLeanEvent<2>,2> * c = b->getBoxAtCoord(coords);
    TS_ASSERT_EQUALS(c, b->getChild(11));
  }


  //-------------------------------------------------------------------------------------
  /** 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 IMDBox<MDLeanEvent<2>,2> ibox_t;

    gbox_t * b = MDEventsTestHelper::makeMDGridBox<2>();
    b->getBoxController()->setSplitThreshold(100);
    b->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
      coord_t centers[2] = {1e-10, 1e-10};
      events.push_back( MDLeanEvent<2>(2.0, 2.0, centers) );
    }
    TS_ASSERT_THROWS_NOTHING( b->addEvents( events ); );


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

    // Dig recursively into the gridded box hierarchies
    std::vector<ibox_t*> boxes;
    size_t expected_depth = 0;
    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);

      // The other boxes have nothing
      TS_ASSERT_EQUALS( boxes[1]->getNPoints(), 0);
      // The other box is a MDBox (it was not split)
      TS_ASSERT( dynamic_cast<box_t*>(boxes[1]) )
    }

    // We went this many levels (and no further) because recursion depth is limited
    TS_ASSERT_EQUALS(boxes[0]->getDepth(), 4);
  }


  //------------------------------------------------------------------------------------------------
  /** This test splits a large number of events, and uses a ThreadPool
   * to use all cores.
   */
  void test_splitAllIfNeeded_usingThreadPool()
  {
    typedef MDGridBox<MDLeanEvent<2>,2> gbox_t;
    typedef MDBox<MDLeanEvent<2>,2> box_t;
    typedef IMDBox<MDLeanEvent<2>,2> ibox_t;

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

    // Make a 1000 events in each sub-box
    size_t num_repeat = 1000;
    if (DODEBUG) num_repeat = 2000;

    Timer tim;
    if (DODEBUG) std::cout << "Adding " << num_repeat*100 << " events...\n";
    MDEventsTestHelper::feedMDBox<2>(b, num_repeat, 10, 0.5, 1.0);
    if (DODEBUG) std::cout << "Adding events done in " << tim.elapsed() << "!\n";

    // Split those boxes in parallel.
    ThreadSchedulerFIFO * ts = new ThreadSchedulerFIFO();
    ThreadPool tp(ts);
    b->splitAllIfNeeded(ts);
    tp.joinAll();

    if (DODEBUG) std::cout << "Splitting events done in " << tim.elapsed() << " sec.\n";

    // Now check the results. Each sub-box should be MDGridBox and have that many events
    std::vector<ibox_t*> boxes = b->getBoxes();
    TS_ASSERT_EQUALS(boxes.size(), 100);
    for (size_t i=0; i<boxes.size(); i++)
    {
      ibox_t * box = boxes[i];
      TS_ASSERT_EQUALS( box->getNPoints(), num_repeat );
      TS_ASSERT( dynamic_cast<gbox_t *>(box) );

      size_t numChildren = box->getNumChildren();
      if (numChildren > 0)
      {
        size_t lastId = box->getChild(0)->getId();
        for (size_t i = 1; i < numChildren; i++)
        {
          TSM_ASSERT_EQUALS("Children IDs need to be sequential!", box->getChild(i)->getId(), lastId+1);
          lastId = box->getChild(i)->getId();
        }
      }

    }

  }




  //------------------------------------------------------------------------------------------------
  /** This test splits a large number of events,
   * for a workspace that is backed by a file (and thus tries to stay below
   * a certain amount of memory used).
   */
  void test_splitAllIfNeeded_fileBacked()
  {
    typedef MDLeanEvent<2> MDE;
    typedef MDGridBox<MDE,2> gbox_t;
    typedef MDBox<MDE,2> box_t;
    typedef IMDBox<MDE,2> ibox_t;

    // Make a fake file-backing for the grid box
    std::string filename = "MDGridBoxTest.nxs";
    ::NeXus::File * file = new ::NeXus::File(filename, NXACC_CREATE);
    file->makeGroup("MDEventWorkspaceTest", "NXentry", 1);
    MDE::prepareNexusData(file, 2000);
    file->close();
    file = new ::NeXus::File(filename, NXACC_RDWR);
    file->openGroup("MDEventWorkspaceTest", "NXentry");
    MDE::openNexusData(file);

    // Create the grid box and make it file-backed.
    gbox_t * b = MDEventsTestHelper::makeMDGridBox<2>();
    BoxController_sptr bc = b->getBoxController();
    bc->setSplitThreshold(100);
    bc->setMaxDepth(4);
    bc->setCacheParameters(1, 1000);
    bc->setFile(file, filename, 0);
    DiskBuffer & dbuf = bc->getDiskBuffer();
    dbuf.setFileLength(0);

    size_t num_repeat = 10;
    if (DODEBUG) num_repeat = 20;
    Timer tim;
    if (DODEBUG) std::cout << "Adding " << num_repeat*10000 << " events...\n";
    MDEventsTestHelper::feedMDBox<2>(b, num_repeat, 100, 0.05, 0.1);
    if (DODEBUG) std::cout << "Adding events done in " << tim.elapsed() << "!\n";

    // Split those boxes in parallel.