Interface for parallel IMDIterators, implemented for MatrixWS

Refs #4782, #4459. Not all builds at this point.

Code/Mantid/Framework/API/inc/MantidAPI/IMDIterator.h

@@ -6,6 +6,7 @@
 //----------------------------------------------------------------------
 #include "MantidAPI/DllConfig.h"
 #include "MantidAPI/IFitFunction.h"
+#include "MantidAPI/IMDWorkspace.h"
 #include "MantidGeometry/MDGeometry/IMDDimension.h"
 #include "MantidGeometry/MDGeometry/MDTypes.h"
 #include "MantidKernel/VMD.h"
@@ -50,19 +51,22 @@ class IMDWorkspace;
   {
   public:
 
+    void setNormalization(Mantid::API::MDNormalization normalization);
+    Mantid::API::MDNormalization getNormalization() const;
+
     /// Get the size of the data (number of entries that will be iterated through)
     virtual size_t getDataSize() const = 0;
 
+    /// Advance to the next cell. If the current cell is the last one in the workspace
+    /// do nothing and return false.
+    virtual bool next() = 0;
+
     /// Is the current position of the iterator valid?
     virtual bool valid() const = 0;
 
     /// Jump to the index^th cell.
     virtual void jumpTo(size_t index) = 0;
 
-    /// Advance to the next cell. If the current cell is the last one in the workspace
-    /// do nothing and return false.
-    virtual bool next() = 0;
-
     /// Advance, skipping a certain number of cells.
     virtual bool next(size_t skip) = 0;
 

Code/Mantid/Framework/API/inc/MantidAPI/IMDWorkspace.h

@@ -12,6 +12,7 @@
 #include <stdarg.h>
 #include "MantidAPI/MDGeometry.h"
 #include "MantidGeometry/MDGeometry/MDImplicitFunction.h"
+#include "MantidAPI/IMDWorkspace.h"
 
 namespace Mantid
 {
@@ -83,20 +84,19 @@ namespace Mantid
       virtual uint64_t getNPoints() const = 0;
 
       /// Creates a new iterator pointing to the first cell in the workspace
-      virtual IMDIterator* createIterator(Mantid::Geometry::MDImplicitFunction * function = NULL) const;
+      virtual std::vector<IMDIterator*> createIterators(size_t suggestedNumCores = 1,
+          Mantid::Geometry::MDImplicitFunction * function = NULL) const = 0;
+
+      /// Creates a new iterator pointing to the first cell in the workspace
+      IMDIterator* createIterator(Mantid::Geometry::MDImplicitFunction * function = NULL) const;
 
       /// Returns the (normalized) signal at a given coordinates
-      virtual signal_t getSignalAtCoord(const coord_t * coords, const Mantid::API::MDNormalization & normalization) const
-      {
-        UNUSED_ARG(coords);
-        UNUSED_ARG(normalization);
-        throw std::runtime_error("getSignalAtCoord() not implemented.");
-      }
+      virtual signal_t getSignalAtCoord(const coord_t * coords, const Mantid::API::MDNormalization & normalization) const = 0;
 
       //-------------------------------------------------------------------------------------------
       /// Returns the signal (normalized by volume) at a given coordinates
       /// @param coords :: coordinate as a VMD vector
-      signal_t getSignalAtCoord(const Mantid::Kernel::VMD & coords, const Mantid::API::MDNormalization & normalization) const
+      signal_t getSignalAtVMD(const Mantid::Kernel::VMD & coords, const Mantid::API::MDNormalization & normalization) const
       {
         return this->getSignalAtCoord(coords.getBareArray(), normalization);
       }

Code/Mantid/Framework/API/inc/MantidAPI/MatrixWorkspace.h

@@ -288,7 +288,8 @@ namespace Mantid
       void saveSpectraMapNexus(::NeXus::File * file, const std::vector<int>& spec,
           const ::NeXus::NXcompression compression = ::NeXus::LZW) const;
 
-      IMDIterator* createIterator(Mantid::Geometry::MDImplicitFunction * function) const;
+      virtual std::vector<IMDIterator*> createIterators(size_t suggestedNumCores = 1,
+          Mantid::Geometry::MDImplicitFunction * function = NULL) const;
 
       void saveInstrumentNexus(::NeXus::File * file) const;
       void loadInstrumentNexus(::NeXus::File * file);

Code/Mantid/Framework/API/inc/MantidAPI/MatrixWorkspaceMDIterator.h

@@ -40,7 +40,8 @@ namespace API
   class DLLExport MatrixWorkspaceMDIterator : public IMDIterator
   {
   public:
-    MatrixWorkspaceMDIterator(const MatrixWorkspace * workspace, Mantid::Geometry::MDImplicitFunction * function);
+    MatrixWorkspaceMDIterator(const MatrixWorkspace * workspace, Mantid::Geometry::MDImplicitFunction * function,
+        size_t beginWI = 0, size_t endWI = size_t(-1));
     virtual ~MatrixWorkspaceMDIterator();
     virtual size_t getDataSize() const;
 
@@ -118,8 +119,12 @@ namespace API
 
     /// Blocksize of workspace
     size_t m_blockSize;
-    /// Number of histograms in workspace
-    size_t m_numHistos;
+
+    /// Workspace index at which the iterator begins
+    size_t m_beginWI;
+
+    /// Workspace index at which the iterator ends
+    size_t m_endWI;
   };
 
 

Code/Mantid/Framework/API/src/IMDWorkspace.cpp

@@ -27,12 +27,17 @@ namespace Mantid
       {
       }
 
-      /**
-        * Default implementation throws NotImplementedError exception.
-        */
-      IMDIterator* IMDWorkspace::createIterator(Mantid::Geometry::MDImplicitFunction * ) const
+      /** Creates a single iterator and returns it.
+       *
+       * This calls createIterators(), a pure virtual method on IMDWorkspace which
+       * has custom implementations for other workspaces.
+       *
+       * @param function :: Implicit function limiting space to look at
+       * @return a single IMDIterator pointer
+       */
+      IMDIterator* IMDWorkspace::createIterator(Mantid::Geometry::MDImplicitFunction * function) const
       {
-        throw Kernel::Exception::NotImplementedError("Iterator is not implemented for this workspace");
+        return this->createIterators(1, function)[0];
       }
   }
 }

Code/Mantid/Framework/API/src/MatrixWorkspace.cpp

@@ -1485,14 +1485,31 @@ namespace Mantid
 
 
     //--------------------------------------------------------------------------------------------
-    /** Create a IMDIterator from this 2D workspace
+    /** Create IMDIterators from this 2D workspace
      *
-     * @param function :: implicit function
-     * @return MatrixWorkspaceMDIterator
+     * @param suggestedNumCores :: split the iterators into this many cores (if threadsafe)
+     * @param function :: implicit function to limit range
+     * @return MatrixWorkspaceMDIterator vector
      */
-    IMDIterator* MatrixWorkspace::createIterator(Mantid::Geometry::MDImplicitFunction * function) const
+    std::vector<IMDIterator*> MatrixWorkspace::createIterators(size_t suggestedNumCores,
+        Mantid::Geometry::MDImplicitFunction * function) const
     {
-      return new MatrixWorkspaceMDIterator(this, function);
+      size_t numCores = suggestedNumCores;
+      if (numCores < 1) numCores = 1;
+      if (!this->threadSafe()) numCores = 1;
+      if (numCores > this->getNumberHistograms())  numCores = this->getNumberHistograms();
+
+      // Create one iterator per core, splitting evenly amongst spectra
+      std::vector<IMDIterator*> out;
+      for (size_t i=0; i<numCores; i++)
+      {
+        size_t beginWI = (i * this->getNumberHistograms()) / numCores;
+        size_t endWI = ((i+1) * this->getNumberHistograms()) / numCores;
+        if (endWI > this->getNumberHistograms())
+          endWI = this->getNumberHistograms();
+        out.push_back(new MatrixWorkspaceMDIterator(this, function, beginWI, endWI));
+      }
+      return out;
     }
 
 

Code/Mantid/Framework/API/src/MatrixWorkspaceMDIterator.cpp

@@ -14,21 +14,40 @@ namespace API
 
   //----------------------------------------------------------------------------------------------
   /** Constructor
+   *
+   * @param workspace :: iterate through this workspace
+   * @param function :: limiting implicit function
+   * @param beginWI :: first workspace index to iterate
+   * @param endWI :: end when you reach this workspace index
    */
-  MatrixWorkspaceMDIterator::MatrixWorkspaceMDIterator(const MatrixWorkspace * workspace, Mantid::Geometry::MDImplicitFunction * function)
+  MatrixWorkspaceMDIterator::MatrixWorkspaceMDIterator(const MatrixWorkspace * workspace,
+      Mantid::Geometry::MDImplicitFunction * function,
+      size_t beginWI, size_t endWI)
   : m_ws(workspace), m_pos(0), m_max(0), m_function(function), m_errorIsCached(false)
   {
     if (!m_ws)
       throw std::runtime_error("MatrixWorkspaceMDIterator::ctor() NULL MatrixWorkspace");
-    m_max = m_ws->size();
     m_center = VMD(2);
     m_isBinnedData = m_ws->isHistogramData();
     m_dimY = m_ws->getDimension(1);
     m_blockSize = m_ws->blocksize();
-    m_numHistos = m_ws->getNumberHistograms();
-    m_workspaceIndex = -1;
+
+    m_beginWI = beginWI;
+    if (m_beginWI >= m_ws->getNumberHistograms())
+      throw std::runtime_error("MatrixWorkspaceMDIterator: Beginning workspace index passed is too high.");
+
+    // End point (handle default)
+    m_endWI = endWI;
+    if (m_endWI > m_ws->getNumberHistograms())
+      m_endWI = m_ws->getNumberHistograms();
+    if (m_endWI < m_beginWI)
+      throw std::runtime_error("MatrixWorkspaceMDIterator: End point is before the start point.");
+
+    m_max = (m_endWI - m_beginWI) * m_blockSize;
     m_xIndex = 0;
-    calcWorkspacePos(0);
+    // Trigger the calculation for the first index
+    m_workspaceIndex = -1;
+    calcWorkspacePos(m_beginWI);
   }
 
   //----------------------------------------------------------------------------------------------
@@ -56,14 +75,15 @@ namespace API
   {
     m_pos = uint64_t(index);
     m_xIndex = m_pos % m_blockSize;
-    size_t newWI = m_pos / m_blockSize;
+    size_t newWI = m_beginWI + (m_pos / m_blockSize);
     calcWorkspacePos(newWI);
   }
 
+  //----------------------------------------------------------------------------------------------
   /// Calculate the workspace index/x index for this iterator position
   inline void MatrixWorkspaceMDIterator::calcWorkspacePos(size_t newWI)
   {
-    if (newWI >= m_numHistos)
+    if (newWI >= m_endWI)
       return;
 
     if (newWI != m_workspaceIndex)

Code/Mantid/Framework/API/test/MatrixWorkspaceMDIteratorTest.h

@@ -17,7 +17,7 @@ class MatrixWorkspaceMDIteratorTest : public CxxTest::TestSuite
 {
 public:
 
-  void test_iterating()
+  boost::shared_ptr<MatrixWorkspace> makeFakeWS()
   {
     boost::shared_ptr<MatrixWorkspace> ws(new WorkspaceTester());
     // Matrix with 4 spectra, 5 bins each
@@ -32,6 +32,12 @@ class MatrixWorkspaceMDIteratorTest : public CxxTest::TestSuite
           ws->dataE(wi)[x] = double((wi*10 + x)*2);
         }
       }
+    return ws;
+  }
+
+  void test_iterating()
+  {
+    boost::shared_ptr<MatrixWorkspace> ws = makeFakeWS();
     IMDIterator * it;
     TS_ASSERT_THROWS_NOTHING( it = ws->createIterator(NULL) );
     TS_ASSERT_EQUALS( it->getDataSize(), 20);
@@ -51,7 +57,34 @@ class MatrixWorkspaceMDIteratorTest : public CxxTest::TestSuite
     TS_ASSERT_DELTA( it->getError(), 22.0, 1e-5);
     TS_ASSERT_DELTA( it->getCenter()[0], 1.5, 1e-5);
     TS_ASSERT_DELTA( it->getCenter()[1], 1.0, 1e-5);
+  }
+
+
+  /** Create a set of iterators that can be applied in parallel */
+  void test_parallel_iterators()
+  {
+    boost::shared_ptr<MatrixWorkspace> ws = makeFakeWS();
+    // The number of output cannot be larger than the number of histograms
+    TS_ASSERT_EQUALS( ws->createIterators(10, NULL).size(), 4);
+
+    // Split in 4 iterators
+    std::vector<IMDIterator*> iterators = ws->createIterators(4, NULL);
+    TS_ASSERT_EQUALS( iterators.size(), 4 );
 
+    for (size_t i=0; i<iterators.size(); i++)
+    {
+      IMDIterator * it = iterators[i];
+      // Only 5 elements per each iterator
+      TS_ASSERT_EQUALS( it->getDataSize(), 5);
+      TS_ASSERT_DELTA( it->getSignal(), double(i)*10 + 0.0, 1e-5);
+      it->next();
+      TS_ASSERT_DELTA( it->getSignal(), double(i)*10 + 1.0, 1e-5);
+      TS_ASSERT_DELTA( it->getError(), double(i)*20 + 2.0, 1e-5);
+      TS_ASSERT( it->next() );
+      TS_ASSERT( it->next() );
+      TS_ASSERT( it->next() );
+      TS_ASSERT( !it->next() );
+    }
   }
 
 

Code/Mantid/Framework/API/test/MatrixWorkspaceTest.h

@@ -602,9 +602,9 @@ class MatrixWorkspaceTest : public CxxTest::TestSuite
         }
       }
     coord_t coords[2] = {0.5, 1.0};
-    TS_ASSERT_DELTA(ws->getSignalAtCoord(coords), 10.0, 1e-5);
+    TS_ASSERT_DELTA(ws->getSignalAtCoord(coords, Mantid::API::NoNormalization), 10.0, 1e-5);
     coords[0] = 1.5;
-    TS_ASSERT_DELTA(ws->getSignalAtCoord(coords), 11.0, 1e-5);
+    TS_ASSERT_DELTA(ws->getSignalAtCoord(coords, Mantid::API::NoNormalization), 11.0, 1e-5);
   }
 
 

Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDEventWorkspace.h

@@ -55,7 +55,8 @@ namespace MDEvents
     virtual uint64_t getNPoints() const;
 
     /// Creates a new iterator pointing to the first cell (box) in the workspace
-    virtual Mantid::API::IMDIterator* createIterator(Mantid::Geometry::MDImplicitFunction * function = NULL) const;
+    virtual std::vector<Mantid::API::IMDIterator*> createIterators(size_t suggestedNumCores = 1,
+        Mantid::Geometry::MDImplicitFunction * function = NULL) const;
 
     /// Returns the (normalized) signal at a given coordinates
     virtual signal_t getSignalAtCoord(const coord_t * coords, const Mantid::API::MDNormalization & normalization) const;

Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspace.h

@@ -64,7 +64,8 @@ namespace MDEvents
       return m_length;
     }
 
-    virtual Mantid::API::IMDIterator* createIterator(Mantid::Geometry::MDImplicitFunction * function = NULL) const;
+    virtual std::vector<Mantid::API::IMDIterator*> createIterators(size_t suggestedNumCores = 1,
+        Mantid::Geometry::MDImplicitFunction * function = NULL) const;
 
     virtual void getLinePlot(const Mantid::Kernel::VMD & start, const Mantid::Kernel::VMD & end,
         Mantid::API::MDNormalization normalize, std::vector<coord_t> & x, std::vector<signal_t> & y, std::vector<signal_t> & e) const;

Code/Mantid/Framework/MDEvents/src/MDEventWorkspace.cpp

@@ -199,16 +199,24 @@ namespace MDEvents
 
   //-----------------------------------------------------------------------------------------------
   /** Creates a new iterator pointing to the first cell (box) in the workspace
+   *
+   * @param suggestedNumCores :: split iterator over this many cores.
    * @param function :: Optional MDImplicitFunction limiting the iterator
+   * @param normalization :: how signal will be normalized
    */
   TMDE(
-  Mantid::API::IMDIterator*  MDEventWorkspace)::createIterator(Mantid::Geometry::MDImplicitFunction * function) const
+  std::vector<Mantid::API::IMDIterator*>  MDEventWorkspace)::createIterators(size_t suggestedNumCores,
+      Mantid::Geometry::MDImplicitFunction * function) const
   {
     // TODO: Should this be leaf only? Depends on most common use case
-    return new MDBoxIterator<MDE,nd>(data, 10000, true, function);
+    // TODO: handle parallelization
+    std::vector<Mantid::API::IMDIterator*> out;
+    out.push_back( new MDBoxIterator<MDE,nd>(data, 10000, true, function) );
+    return out;
   }
 
 
+
   //-----------------------------------------------------------------------------------------------
   /** Returns the (normalized) signal at a given coordinates
    *

Code/Mantid/Framework/MDEvents/src/MDHistoWorkspace.cpp

@@ -351,9 +351,17 @@ namespace MDEvents
 
 
   //----------------------------------------------------------------------------------------------
-  /// Creates a new iterator pointing to the first cell in the workspace
-  Mantid::API::IMDIterator* MDHistoWorkspace::createIterator(Mantid::Geometry::MDImplicitFunction * function) const
+  /** Create IMDIterators from this MDHistoWorkspace
+   *
+   * @param suggestedNumCores :: split the iterators into this many cores (if threadsafe)
+   * @param function :: implicit function to limit range
+   * @return MatrixWorkspaceMDIterator vector
+   */
+
+  std::vector<Mantid::API::IMDIterator*> MDHistoWorkspace::createIterators(size_t suggestedNumCores,
+      Mantid::Geometry::MDImplicitFunction * function) const
   {
+    //TODO:
     return new MDHistoWorkspaceIterator(this,function);
   }