Refs #4036 all standard binary operations and boolean ops on MDHistoWorkspace

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

@@ -75,6 +75,12 @@ namespace MDEvents
     MDHistoWorkspace & operator/=(const MDHistoWorkspace & b);
     MDHistoWorkspace & operator/=(const Mantid::DataObjects::WorkspaceSingleValue & b);
 
+    MDHistoWorkspace & operator&=(const MDHistoWorkspace & b);
+    MDHistoWorkspace & operator|=(const MDHistoWorkspace & b);
+    MDHistoWorkspace & operator^=(const MDHistoWorkspace & b);
+
+    void operatorNot();
+
 
 
     // --------------------------------------------------------------------------------------------

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

@@ -324,14 +324,22 @@ namespace MDEvents
     return out;
   }
 
+  //----------------------------------------------------------------------------------------------
+  /** Check if the two workspace's sizes match (for comparison or element-by-element operation
+   *
+   * @param other :: the workspace to compare to
+   * @param operation :: descriptive string (for the error message)
+   * @throw an error if they don't match
+   */
   void MDHistoWorkspace::checkWorkspaceSize(const MDHistoWorkspace & other, std::string operation)
   {
     if (other.getNumDims() != this->getNumDims())
-      throw std::invalid_argument("Cannot perform the " + operation + " operation on this workspace. The number of dimensions does not match.");
+      throw std::invalid_argument("Cannot perform the " + operation + " operation on this MDHistoWorkspace. The number of dimensions does not match.");
     if (other.m_length != this->m_length)
-      throw std::invalid_argument("Cannot perform the " + operation + " operation on this workspace. The length of the signals vector does not match.");
+      throw std::invalid_argument("Cannot perform the " + operation + " operation on this MDHistoWorkspace. The length of the signals vector does not match.");
   }
 
+
   //----------------------------------------------------------------------------------------------
   /** Perform the += operation, element-by-element, for two MDHistoWorkspace's
    *
@@ -366,6 +374,220 @@ namespace MDEvents
     return *this;
   }
 
+  //----------------------------------------------------------------------------------------------
+  /** Perform the -= operation, element-by-element, for two MDHistoWorkspace's
+   *
+   * @param b :: workspace on the RHS of the operation
+   * @return *this after operation */
+  MDHistoWorkspace & MDHistoWorkspace::operator-=(const MDHistoWorkspace & b)
+  {
+    checkWorkspaceSize(b, "-=");
+    for (size_t i=0; i<m_length; ++i)
+    {
+      m_signals[i] -= b.m_signals[i];
+      m_errorsSquared[i] += b.m_errorsSquared[i];
+    }
+    return *this;
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** Perform the -= operation with a scalar as the RHS argument
+   *
+   * @param b :: WorkspaceSingleValue (signal and error) as the RHS argument
+   * @return *this after operation */
+  MDHistoWorkspace & MDHistoWorkspace::operator-=(const Mantid::DataObjects::WorkspaceSingleValue & b)
+  {
+    signal_t signal = b.dataY(0)[0];
+    signal_t errorSquared = b.dataE(0)[0];
+    errorSquared *= errorSquared;
+    for (size_t i=0; i<m_length; ++i)
+    {
+      m_signals[i] -= signal;
+      m_errorsSquared[i] += errorSquared;
+    }
+    return *this;
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** Perform the *= operation, element-by-element, for two MDHistoWorkspace's
+   *
+   * Error propagation of \f$ f = a * b \f$  is given by:
+   * \f$ df^2 = f^2 * (da^2 / a^2 + db^2 / b^2) \f$
+   *
+   * @param b_ws :: workspace on the RHS of the operation
+   * @return *this after operation */
+  MDHistoWorkspace & MDHistoWorkspace::operator*=(const MDHistoWorkspace & b_ws)
+  {
+    checkWorkspaceSize(b_ws, "*=");
+    for (size_t i=0; i<m_length; ++i)
+    {
+      signal_t a = m_signals[i];
+      signal_t da2 = m_errorsSquared[i];
+
+      signal_t b = b_ws.m_signals[i];
+      signal_t db2 = b_ws.m_errorsSquared[i];
+
+      signal_t f = a * b;
+      signal_t df2 = (f * f) * ( da2 / (a*a) + db2 / (b*b) );
+
+      m_signals[i] = f;
+      m_errorsSquared[i] = df2;
+    }
+    return *this;
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** Perform the *= operation with a scalar as the RHS argument
+   *
+   * Error propagation of \f$ f = a * b \f$  is given by:
+   * \f$ df^2 = f^2 * (da^2 / a^2 + db^2 / b^2) \f$
+   *
+   * @param b_ws :: WorkspaceSingleValue (signal and error) as the RHS argument
+   * @return *this after operation */
+  MDHistoWorkspace & MDHistoWorkspace::operator*=(const Mantid::DataObjects::WorkspaceSingleValue & b_ws)
+  {
+    signal_t b = b_ws.dataY(0)[0];
+    signal_t db2 = b_ws.dataE(0)[0] * b_ws.dataE(0)[0];
+    signal_t db2_relative = db2 / (b*b);
+    for (size_t i=0; i<m_length; ++i)
+    {
+      signal_t a = m_signals[i];
+      signal_t da2 = m_errorsSquared[i];
+
+      signal_t f = a * b;
+      signal_t df2 = (f * f) * ( da2 / (a*a) + db2_relative );
+
+      m_signals[i] = f;
+      m_errorsSquared[i] = df2;
+    }
+    return *this;
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** Perform the /= operation, element-by-element, for two MDHistoWorkspace's
+   *
+   * Error propagation of \f$ f = a / b \f$  is given by:
+   * \f$ df^2 = f^2 * (da^2 / a^2 + db^2 / b^2) \f$
+   *
+   * @param b_ws :: workspace on the RHS of the operation
+   * @return *this after operation */
+  MDHistoWorkspace & MDHistoWorkspace::operator/=(const MDHistoWorkspace & b_ws)
+  {
+    checkWorkspaceSize(b_ws, "/=");
+    for (size_t i=0; i<m_length; ++i)
+    {
+      signal_t a = m_signals[i];
+      signal_t da2 = m_errorsSquared[i];
+
+      signal_t b = b_ws.m_signals[i];
+      signal_t db2 = b_ws.m_errorsSquared[i];
+
+      signal_t f = a / b;
+      signal_t df2 = (f * f) * ( da2 / (a*a) + db2 / (b*b) );
+
+      m_signals[i] = f;
+      m_errorsSquared[i] = df2;
+    }
+    return *this;
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** Perform the /= operation with a scalar as the RHS argument
+   *
+   * Error propagation of \f$ f = a / b \f$  is given by:
+   * \f$ df^2 = f^2 * (da^2 / a^2 + db^2 / b^2) \f$
+   *
+   * @param b_ws :: WorkspaceSingleValue (signal and error) as the RHS argument
+   * @return *this after operation */
+  MDHistoWorkspace & MDHistoWorkspace::operator/=(const Mantid::DataObjects::WorkspaceSingleValue & b_ws)
+  {
+    signal_t b = b_ws.dataY(0)[0];
+    signal_t db2 = b_ws.dataE(0)[0] * b_ws.dataE(0)[0];
+    signal_t db2_relative = db2 / (b*b);
+    for (size_t i=0; i<m_length; ++i)
+    {
+      signal_t a = m_signals[i];
+      signal_t da2 = m_errorsSquared[i];
+
+      signal_t f = a / b;
+      signal_t df2 = (f * f) * ( da2 / (a*a) + db2_relative );
+
+      m_signals[i] = f;
+      m_errorsSquared[i] = df2;
+    }
+    return *this;
+  }
+
+
+  //----------------------------------------------------------------------------------------------
+  /** A boolean &= (and) operation, element-by-element, for two MDHistoWorkspace's.
+   *
+   * 0.0 is "false", all other values are "true". All errors are set to 0.
+   *
+   * @param b :: workspace on the RHS of the operation
+   * @return *this after operation */
+  MDHistoWorkspace & MDHistoWorkspace::operator&=(const MDHistoWorkspace & b)
+  {
+    checkWorkspaceSize(b, "&= (and)");
+    for (size_t i=0; i<m_length; ++i)
+    {
+      m_signals[i] = m_signals[i] && b.m_signals[i];
+      m_errorsSquared[i] = 0;
+    }
+    return *this;
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** A boolean |= (or) operation, element-by-element, for two MDHistoWorkspace's.
+   *
+   * 0.0 is "false", all other values are "true". All errors are set to 0.
+   *
+   * @param b :: workspace on the RHS of the operation
+   * @return *this after operation */
+  MDHistoWorkspace & MDHistoWorkspace::operator|=(const MDHistoWorkspace & b)
+  {
+    checkWorkspaceSize(b, "|= (or)");
+    for (size_t i=0; i<m_length; ++i)
+    {
+      m_signals[i] = m_signals[i] || b.m_signals[i];
+      m_errorsSquared[i] = 0;
+    }
+    return *this;
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** A boolean ^= (xor) operation, element-by-element, for two MDHistoWorkspace's.
+   *
+   * 0.0 is "false", all other values are "true". All errors are set to 0.
+   *
+   * @param b :: workspace on the RHS of the operation
+   * @return *this after operation */
+  MDHistoWorkspace & MDHistoWorkspace::operator^=(const MDHistoWorkspace & b)
+  {
+    checkWorkspaceSize(b, "^= (xor)");
+    for (size_t i=0; i<m_length; ++i)
+    {
+      m_signals[i] = bool(m_signals[i]) ^ bool(b.m_signals[i]);
+      m_errorsSquared[i] = 0;
+    }
+    return *this;
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** A boolean not operation, performed in-place.
+   * All errors are set to 0.
+   *
+   * 0.0 is "false", all other values are "true". All errors are set to 0.
+   */
+  void MDHistoWorkspace::operatorNot()
+  {
+    for (size_t i=0; i<m_length; ++i)
+    {
+      m_signals[i] = !bool(m_signals[i]);
+      m_errorsSquared[i] = 0;
+    }
+  }
+
 
 } // namespace Mantid
 } // namespace MDEvents

Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceTest.h

@@ -14,6 +14,7 @@
 #include <iomanip>
 #include <iostream>
 #include "MantidAPI/IMDWorkspace.h"
+#include "MantidDataObjects/WorkspaceSingleValue.h"
 
 using namespace Mantid::MDEvents;
 using namespace Mantid::Geometry;
@@ -32,16 +33,25 @@ public:
   MDHistoWorkspace_sptr two;
   MDHistoWorkspace_sptr three;
 
-  MDHistoWorkspace_sptr ma
-
+  /** Create some fake workspace, 5x5, with 2.0 and 3.0 each */
   MDHistoWorkspaceTest()
   {
-    MDEventsTestHelper::makeFakeMDHistoWorkspace(
-    two.
+    two = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.0);
+    three = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 3.0);
   }
 
+  /** Check that a workspace has the right signal/error*/
+  void checkWorkspace(MDHistoWorkspace_sptr ws, double expectedSignal, double expectedErrorSquared)
+  {
+    for (size_t i=0; i < ws->getNPoints(); i++)
+    {
+      TS_ASSERT_DELTA( ws->getSignalAt(i), expectedSignal, 1e-5 );
+      TS_ASSERT_DELTA( ws->getErrorAt(i), expectedErrorSquared, 1e-5 );
+    }
+  }
 
 
+  //--------------------------------------------------------------------------------------
   void test_constructor()
   {
     MDHistoDimension_sptr dimX(new MDHistoDimension("X", "x", "m", -10, 10, 5));
@@ -364,7 +374,7 @@ public:
   void test_getSignalAtCoord()
   {
     // 2D workspace with signal[i] = i (linear index)
-    MDHistoWorkspace_sptr ws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0, 1.0, 2, 10);
+    MDHistoWorkspace_sptr ws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0, 2, 10);
     for (size_t i=0; i<100; i++)
       ws->setSignalAt(i, double(i));
     IMDWorkspace_sptr iws(ws);
@@ -379,6 +389,133 @@ public:
     TS_ASSERT( boost::math::isnan(iws->getSignalAtCoord(VMD(3.5, 10.02)) ) );
   }
 
+  //--------------------------------------------------------------------------------------
+  void test_plus_ws()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.5 /*errorSquared*/);
+    MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 3.5 /*errorSquared*/);
+    *a += *b;
+    checkWorkspace(a, 5.0, 6.0);
+  }
+
+  void test_plus_scalar()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.5 /*errorSquared*/);
+    WorkspaceSingleValue b(3.0, sqrt(3.5) );
+    *a += b;
+    checkWorkspace(a, 5.0, 6.0);
+  }
+
+  //--------------------------------------------------------------------------------------
+  void test_minus_ws()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 2.5 /*errorSquared*/);
+    MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 3.5 /*errorSquared*/);
+    *a -= *b;
+    checkWorkspace(a, 1.0, 6.0);
+  }
+
+  void test_minus_scalar()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 2.5 /*errorSquared*/);
+    WorkspaceSingleValue b(2.0, sqrt(3.5) );
+    *a -= b;
+    checkWorkspace(a, 1.0, 6.0);
+  }
+
+  //--------------------------------------------------------------------------------------
+  void test_times_ws()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.0 /*errorSquared*/);
+    MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 3.0 /*errorSquared*/);
+    *a *= *b;
+    checkWorkspace(a, 6.0, 36. * (.5 + 1./3.));
+  }
+
+  //--------------------------------------------------------------------------------------
+  void test_times_scalar()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.0 /*errorSquared*/);
+    WorkspaceSingleValue b(3.0, sqrt(3.0) );
+    *a *= b;
+    checkWorkspace(a, 6.0, 36. * (.5 + 1./3.));
+    // Scalar without error
+    MDHistoWorkspace_sptr d = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.0 /*errorSquared*/);
+    WorkspaceSingleValue e(3.0, 0);
+    *d *= e;
+    checkWorkspace(d, 6.0, 9 * 2.0);
+  }
+
+  //--------------------------------------------------------------------------------------
+  void test_divide_ws()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 3.0 /*errorSquared*/);
+    MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.0 /*errorSquared*/);
+    *a /= *b;
+    checkWorkspace(a, 1.5, 1.5 * 1.5 * (.5 + 1./3.));
+  }
+
+  //--------------------------------------------------------------------------------------
+  void test_divide_scalar()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 3.0 /*errorSquared*/);
+    WorkspaceSingleValue b(2.0, sqrt(2.0) );
+    *a /= b;
+    checkWorkspace(a, 1.5, 1.5 * 1.5 * (.5 + 1./3.));
+  }
+
+  //--------------------------------------------------------------------------------------
+  void test_boolean_and()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.23, 2, 5, 10.0, 3.0);
+    MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.34, 2, 5, 10.0, 2.0);
+    MDHistoWorkspace_sptr c = MDEventsTestHelper::makeFakeMDHistoWorkspace(0.00, 2, 5, 10.0, 2.0);
+    *a &= *b;
+    checkWorkspace(a, 1.0, 0.0);
+    *b &= *c;
+    checkWorkspace(b, 0.0, 0.0);
+  }
+
+  //--------------------------------------------------------------------------------------
+  void test_boolean_or()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.23, 2, 5, 10.0, 3.0);
+    MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.34, 2, 5, 10.0, 2.0);
+    MDHistoWorkspace_sptr c = MDEventsTestHelper::makeFakeMDHistoWorkspace(0.00, 2, 5, 10.0, 2.0);
+    *a |= *b;
+    checkWorkspace(a, 1.0, 0.0);
+    *b |= *c;
+    checkWorkspace(b, 1.0, 0.0);
+    *c |= *c;
+    checkWorkspace(c, 0.0, 0.0);
+  }
+
+  //--------------------------------------------------------------------------------------
+  void test_boolean_xor()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.23, 2, 5, 10.0, 3.0);
+    MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.34, 2, 5, 10.0, 2.0);
+    MDHistoWorkspace_sptr c = MDEventsTestHelper::makeFakeMDHistoWorkspace(0.00, 2, 5, 10.0, 2.0);
+    *a ^= *b;
+    checkWorkspace(a, 0.0, 0.0);
+    *b ^= *c;
+    checkWorkspace(b, 1.0, 0.0);
+    *c ^= *c;
+    checkWorkspace(c, 0.0, 0.0);
+  }
+
+  //--------------------------------------------------------------------------------------
+  void test_boolean_operatorNot()
+  {
+    MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.23, 2, 5, 10.0, 3.0);
+    MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(0.00, 2, 5, 10.0, 2.0);
+    a->operatorNot();
+    checkWorkspace(a, 0.0, 0.0);
+    b->operatorNot();
+    checkWorkspace(b, 1.0, 0.0);
+  }
+
+
 };
 
 

Code/Mantid/Framework/TestHelpers/inc/MantidTestHelpers/MDEventsTestHelper.h

@@ -43,7 +43,8 @@ namespace MDEventsTestHelper
   DLL_TESTHELPERS Mantid::MDEvents::MDEventWorkspace3Lean::sptr makeFileBackedMDEW(std::string wsName, bool fileBacked);
 
   /// Make a fake n-dimensional MDHistoWorkspace
-  DLL_TESTHELPERS Mantid::MDEvents::MDHistoWorkspace_sptr makeFakeMDHistoWorkspace(double signal, size_t numDims, size_t numBins = 10, double max = 10.0);
+  DLL_TESTHELPERS Mantid::MDEvents::MDHistoWorkspace_sptr makeFakeMDHistoWorkspace(double signal, size_t numDims, size_t numBins = 10, double max = 10.0,
+      double errorSquared=1.0);
 
 
   //-------------------------------------------------------------------------------------

Code/Mantid/Framework/TestHelpers/src/MDEventsTestHelper.cpp

@@ -205,7 +205,7 @@ namespace MDEventsTestHelper
    * @return the MDHisto
    */
   Mantid::MDEvents::MDHistoWorkspace_sptr makeFakeMDHistoWorkspace(double signal, size_t numDims, size_t numBins,
-      double max)
+      double max, double errorSquared)
   {
     Mantid::MDEvents::MDHistoWorkspace * ws = NULL;
     if (numDims ==1)
@@ -236,7 +236,7 @@ namespace MDEventsTestHelper
           );
     }
     Mantid::MDEvents::MDHistoWorkspace_sptr ws_sptr(ws);
-    ws_sptr->setTo(signal, signal);
+    ws_sptr->setTo(signal, errorSquared);
     return ws_sptr;
   }