From 0c6ba424ec177e2b717093e67a3b7a88f95bdd1a Mon Sep 17 00:00:00 2001
From: Janik Zikovsky <zikovskyjl@ornl.gov>
Date: Thu, 17 Mar 2011 18:37:39 +0000
Subject: [PATCH] Refs #2012: Double speed of gridding MDEvents by
 parallelizing the MDBox splitting portion.

---
 .../MDEvents/inc/MantidMDEvents/MDGridBox.h   |  5 +-
 .../Framework/MDEvents/src/MDGridBox.cpp      | 60 +++++++++++++-----
 .../Framework/MDEvents/test/MDGridBoxTest.h   | 62 ++++++++++++++++---
 3 files changed, 101 insertions(+), 26 deletions(-)

diff --git a/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDGridBox.h b/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDGridBox.h
index 728e83b4bfc..0074e7094f4 100644
--- a/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDGridBox.h
+++ b/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDGridBox.h
@@ -9,6 +9,7 @@
 #include "MantidMDEvents/MDBox.h"
 #include "MantidMDEvents/MDDimensionExtents.h"
 #include "MantidMDEvents/MDEvent.h"
+#include "MantidKernel/ThreadScheduler.h"
 
 namespace Mantid
 {
@@ -53,9 +54,9 @@ namespace MDEvents
 
     size_t addManyEvents(const std::vector<MDE> & events, Mantid::API::Progress * prog = NULL);
 
-    void splitContents(size_t index);
+    void splitContents(size_t index, Kernel::ThreadScheduler * ts = NULL);
 
-    void splitAllIfNeeded();
+    void splitAllIfNeeded(Kernel::ThreadScheduler * ts = NULL);
 
     // ======================= Testing/Debugging Methods =================
     /** For testing: get 9a copy of) the vector of boxes */
diff --git a/Code/Mantid/Framework/MDEvents/src/MDGridBox.cpp b/Code/Mantid/Framework/MDEvents/src/MDGridBox.cpp
index 5c1b0284f38..993361e9b79 100644
--- a/Code/Mantid/Framework/MDEvents/src/MDGridBox.cpp
+++ b/Code/Mantid/Framework/MDEvents/src/MDGridBox.cpp
@@ -1,6 +1,7 @@
 #include "MantidAPI/Progress.h"
 #include "MantidKernel/FunctionTask.h"
 #include "MantidKernel/Task.h"
+#include "MantidKernel/FunctionTask.h"
 #include "MantidKernel/Timer.h"
 #include "MantidKernel/ThreadPool.h"
 #include "MantidKernel/ThreadScheduler.h"
@@ -224,9 +225,11 @@ namespace MDEvents
    *
    * @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)
+  void MDGridBox)::splitContents(size_t index, ThreadScheduler * ts)
   {
     // You can only split it if it is a MDBox (not MDGridBox).
     MDBox<MDE, nd> * box = dynamic_cast<MDBox<MDE, nd> *>(boxes[index]);
@@ -237,17 +240,23 @@ namespace MDEvents
     delete boxes[index];
     // And now we have a gridded box instead of a boring old regular box.
     boxes[index] = gridbox;
+
+    if (ts)
+    {
+      // Create a task to split the newly create MDGridBox.
+      ts->push(new FunctionTask(boost::bind(&MDGridBox<MDE,nd>::splitAllIfNeeded, &*gridbox, ts) ) );
+    }
   }
 
   //-----------------------------------------------------------------------------------------------
   /** Goes through all the sub-boxes and splits them if they contain
    * enough events to be worth it.
    *
-   * @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 to do it serially.
    */
   TMDE(
-  void MDGridBox)::splitAllIfNeeded()
+  void MDGridBox)::splitAllIfNeeded(ThreadScheduler * ts)
   {
     for (size_t i=0; i < numBoxes; ++i)
     {
@@ -257,15 +266,25 @@ namespace MDEvents
         // Plain MD-Box. Does it need to split?
         if (this->m_BoxController->willSplit(box->getNPoints(), box->getDepth() ))
         {
-          //std::cout << "Splitting box at depth " << box->getDepth() << "\n";
           // The MDBox needs to split into a grid box.
-          MDGridBox<MDE, nd> * gridBox = new MDGridBox<MDE, nd>(box);
-          // Replace in the array
-          boxes[i] = gridBox;
-          // Delete the old box
-          delete box;
-          // Now recursively check if this NEW grid box's contents should be split too
-          gridBox->splitAllIfNeeded();
+          if (!ts)
+          {
+            // ------ Perform split serially (no ThreadPool) ------
+            MDGridBox<MDE, nd> * gridBox = new MDGridBox<MDE, nd>(box);
+            // Replace in the array
+            boxes[i] = gridBox;
+            // Delete the old box
+            delete box;
+            // Now recursively check if this NEW grid box's contents should be split too
+            gridBox->splitAllIfNeeded(NULL);
+          }
+          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 FunctionTask(boost::bind(&MDGridBox<MDE,nd>::splitContents, &*this, i, ts) ) );
+          }
         }
       }
       else
@@ -275,7 +294,7 @@ namespace MDEvents
         if (gridBox)
         {
           // Now recursively check if this old grid box's contents should be split too
-          gridBox->splitAllIfNeeded();
+          gridBox->splitAllIfNeeded(ts);
         }
       }
     }
@@ -411,7 +430,7 @@ namespace MDEvents
       //Since the costs are not known ahead of time, use a simple FIFO buffer.
       ThreadScheduler * ts = new ThreadSchedulerFIFO();
       // Create the threadpool
-      ThreadPool tp(ts, 1);
+      ThreadPool tp(ts);
 
       // Do 'numTasksPerBlock' tasks with 'eventsPerTask' events in each one.
       for (size_t i = 0; i < numTasksPerBlock; i++)
@@ -440,11 +459,18 @@ namespace MDEvents
       tp.joinAll();
 //      std::cout << "... block took " << tim.elapsed() << " secs.\n";
 
+
+      //Create a threadpool for splitting.
+      ThreadScheduler * ts_splitter = new ThreadSchedulerFIFO();
+      ThreadPool tp_splitter(ts_splitter);
+
       //Now, shake out all the sub boxes and split those if needed
-//      std::cout << "Starting splitAllIfNeeded().\n";
+      //std::cout << "\nStarting splitAllIfNeeded().\n";
       if (prog) prog->report("Splitting MDBox'es.");
-      this->splitAllIfNeeded();
-//      std::cout << "... splitAllIfNeeded() took " << tim.elapsed() << " secs.\n";
+
+      this->splitAllIfNeeded(ts_splitter);
+      tp_splitter.joinAll();
+      //std::cout << "\n... splitAllIfNeeded() took " << tim.elapsed() << " secs.\n";
     }
 
     // Refresh the counts, now that we are all done.
diff --git a/Code/Mantid/Framework/MDEvents/test/MDGridBoxTest.h b/Code/Mantid/Framework/MDEvents/test/MDGridBoxTest.h
index 84ffcf7c357..4cf0d60a4b3 100644
--- a/Code/Mantid/Framework/MDEvents/test/MDGridBoxTest.h
+++ b/Code/Mantid/Framework/MDEvents/test/MDGridBoxTest.h
@@ -539,8 +539,9 @@ public:
     }
     TS_ASSERT_THROWS_NOTHING( b->addEvents( events ); );
 
+
     // Split into sub-grid boxes
-    TS_ASSERT_THROWS_NOTHING( b->splitAllIfNeeded(); )
+    TS_ASSERT_THROWS_NOTHING( b->splitAllIfNeeded(NULL); )
 
     // Dig recursively into the gridded box hierarchies
     std::vector<ibox_t*> boxes;
@@ -571,13 +572,60 @@ public:
 
     // 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<MDEvent<2>,2> gbox_t;
+    typedef MDBox<MDEvent<2>,2> box_t;
+    typedef IMDBox<MDEvent<2>,2> ibox_t;
+
+    gbox_t * b = makeMDGridBox<2>();
+    b->getBoxController()->m_SplitThreshold = 100;
+    b->getBoxController()->m_maxDepth = 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";
 
+    std::vector< MDEvent<2> > events;
+    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};
+        for (size_t i=0; i < num_repeat; i++)
+        {
+          // Make an event in the middle of each box
+          events.push_back( MDEvent<2>(2.0, 2.0, centers) );
+        }
+      }
+    TS_ASSERT_THROWS_NOTHING( b->addEvents( events ); );
+    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++)
+    {
+      TS_ASSERT_EQUALS( boxes[i]->getNPoints(), num_repeat );
+      TS_ASSERT( dynamic_cast<gbox_t *>(boxes[i]) );
+    }
 
-//    // Get the right totals again
-//    b->refreshCache();
-//    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);
   }
 
   //-------------------------------------------------------------------------------------
@@ -649,7 +697,7 @@ public:
     if (!DODEBUG) return;
 
     ProgressText * prog = new ProgressText(0.0, 1.0, 10, true);
-    prog->setNotifyStep(0.00001); //Notify always
+    prog->setNotifyStep(0.5); //Notify more often
 
     typedef MDGridBox<MDEvent<2>,2> box_t;
     box_t * b = makeMDGridBox<2>();
-- 
GitLab