BinaryOperation.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/BinaryOperation.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/EventList.h"
#include "MantidAPI/FrameworkManager.h"
#include "MantidAPI/MemoryManager.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/WorkspaceOpOverloads.h"
#include "MantidGeometry/IDetector.h"
#include "MantidKernel/Timer.h"
#include "MantidDataObjects/WorkspaceSingleValue.h"

#include <boost/make_shared.hpp>

using namespace Mantid::Geometry;
using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::DataObjects;
using std::size_t;

namespace Mantid {
namespace Algorithms {
BinaryOperation::BinaryOperation()
    : API::Algorithm(), m_lhs(), m_elhs(), m_rhs(), m_erhs(), m_out(), m_eout(),
      m_AllowDifferentNumberSpectra(false), m_ClearRHSWorkspace(false),
      m_matchXSize(false), m_flipSides(false), m_keepEventWorkspace(false),
      m_useHistogramForRhsEventWorkspace(false),
      m_do2D_even_for_SingleColumn_on_rhs(false), m_indicesToMask(),
      m_progress(NULL) {}

BinaryOperation::~BinaryOperation() {
  if (m_progress)
    delete m_progress;
}

/** Initialisation method.
 *  Defines input and output workspaces
 *
 */
void BinaryOperation::init() {
  declareProperty(
      new WorkspaceProperty<MatrixWorkspace>(inputPropName1(), "",
                                             Direction::Input),
      "The name of the input workspace on the left hand side of the operation");
  declareProperty(new WorkspaceProperty<MatrixWorkspace>(inputPropName2(), "",
                                                         Direction::Input),
                  "The name of the input workspace on the right hand side of "
                  "the operation");
  declareProperty(new WorkspaceProperty<MatrixWorkspace>(outputPropName(), "",
                                                         Direction::Output),
                  "The name to call the output workspace");
  declareProperty(
      new PropertyWithValue<bool>("AllowDifferentNumberSpectra", false,
                                  Direction::Input),
      "Are workspaces with different number of spectra allowed? "
      "For example, the LHSWorkspace might have one spectrum per detector, "
      "but the RHSWorkspace could have its spectra averaged per bank. If true, "
      "then matching between the LHS and RHS spectra is performed (all "
      "detectors "
      "in a LHS spectrum have to be in the corresponding RHS) in order to "
      "apply the RHS spectrum to the LHS.");

  declareProperty(
      new PropertyWithValue<bool>("ClearRHSWorkspace", false, Direction::Input),
      "For EventWorkspaces only. This will clear out event lists "
      "from the RHS workspace as the binary operation is applied. "
      "This can prevent excessive memory use, e.g. when subtracting "
      "an EventWorkspace from another: memory use will be approximately "
      "constant instead of increasing by 50%. At completion, the RHS workspace "
      "will be empty.");
}

//--------------------------------------------------------------------------------------------
/** Special handling for 1-WS and 1/WS.
 *
 * @return true if the operation was handled; exec() should then return
 */
bool BinaryOperation::handleSpecialDivideMinus() {
  // Is the LHS operand a single number?
  WorkspaceSingleValue_const_sptr lhs_singleVal =
      boost::dynamic_pointer_cast<const WorkspaceSingleValue>(m_lhs);
  WorkspaceSingleValue_const_sptr rhs_singleVal =
      boost::dynamic_pointer_cast<const WorkspaceSingleValue>(m_rhs);

  if (lhs_singleVal) {
    MatrixWorkspace_sptr out = getProperty("OutputWorkspace");
    if (this->name() == "Divide" && !bool(rhs_singleVal)) {
      // x / workspace = Power(workspace, -1) * x
      // workspace ^ -1
      IAlgorithm_sptr pow = this->createChildAlgorithm("Power", 0.0, 0.5, true);
      pow->setProperty("InputWorkspace",
                       boost::const_pointer_cast<MatrixWorkspace>(m_rhs));
      pow->setProperty("Exponent", -1.0);
      pow->setProperty("OutputWorkspace", out);
      pow->executeAsChildAlg();
      out = pow->getProperty("OutputWorkspace");

      // Multiply by x
      IAlgorithm_sptr mult =
          this->createChildAlgorithm("Multiply", 0.5, 1.0, true);
      mult->setProperty(inputPropName1(), out); //(workspace^-1)
      mult->setProperty(inputPropName2(),
                        boost::const_pointer_cast<MatrixWorkspace>(
                            m_lhs)); // (1.0) or other number
      mult->setProperty(outputPropName(), out);
      mult->executeAsChildAlg();
      out = mult->getProperty("OutputWorkspace");
      setProperty("OutputWorkspace", out);
      return true;
    } else if (this->name() == "Minus") {
      // x - workspace = x + (workspace * -1)
      MatrixWorkspace_sptr minusOne =
          WorkspaceFactory::Instance().create("WorkspaceSingleValue", 1, 1, 1);
      minusOne->dataY(0)[0] = -1.0;
      minusOne->dataE(0)[0] = 0.0;

      // workspace * -1
      IAlgorithm_sptr mult =
          this->createChildAlgorithm("Multiply", 0.0, 0.5, true);
      mult->setProperty(inputPropName1(),
                        boost::const_pointer_cast<MatrixWorkspace>(m_rhs));
      mult->setProperty(inputPropName2(), minusOne);
      mult->setProperty("OutputWorkspace", out);
      mult->executeAsChildAlg();
      out = mult->getProperty("OutputWorkspace");

      // Multiply by x
      IAlgorithm_sptr plus = this->createChildAlgorithm("Plus", 0.5, 1.0, true);
      plus->setProperty(inputPropName1(), out); //(workspace^-1)
      plus->setProperty(inputPropName2(),
                        boost::const_pointer_cast<MatrixWorkspace>(
                            m_lhs)); // (1.0) or other number
      plus->setProperty(outputPropName(), out);
      plus->executeAsChildAlg();
      out = plus->getProperty("OutputWorkspace");
      setProperty("OutputWorkspace", out);
      return true;
    }

  } // lhs_singleVal

  // Process normally
  return false;
}

//--------------------------------------------------------------------------------------------
/** Executes the algorithm. Will call execEvent() if appropriate.
 *
 *  @throw runtime_error Thrown if algorithm cannot execute
 */
void BinaryOperation::exec() {
  // get input workspace, dynamic cast not needed
  m_lhs = getProperty(inputPropName1());
  m_rhs = getProperty(inputPropName2());
  m_AllowDifferentNumberSpectra = getProperty("AllowDifferentNumberSpectra");

  // Special handling for 1-WS and 1/WS.
  if (this->handleSpecialDivideMinus())
    return;

  // Cast to EventWorkspace pointers
  m_elhs = boost::dynamic_pointer_cast<const EventWorkspace>(m_lhs);
  m_erhs = boost::dynamic_pointer_cast<const EventWorkspace>(m_rhs);

  // We can clear the RHS workspace if it is an event,
  //  and we are not doing mismatched spectra (in which case you might clear it
  //  too soon!)
  //  and lhs is not rhs.
  //  and out is not rhs.
  m_ClearRHSWorkspace = getProperty("ClearRHSWorkspace");
  if (m_ClearRHSWorkspace) {
    if (m_AllowDifferentNumberSpectra || (!m_erhs) || (m_rhs == m_lhs) ||
        (m_out == m_rhs)) {
      // std::cout << "m_ClearRHSWorkspace = false\n";
      m_ClearRHSWorkspace = false;
    }
  }

  // Get the output workspace
  m_out = getProperty(outputPropName());
  m_eout = boost::dynamic_pointer_cast<EventWorkspace>(m_out);

  // Make a check of what will be needed to setup the workspaces, based on the
  // input types.
  this->checkRequirements();

  if (m_flipSides) {
    // Flip the workspaces left and right
    std::swap(m_lhs, m_rhs);
    std::swap(m_elhs, m_erhs);
  }

  // Check that the input workspaces are compatible
  if (!checkCompatibility(m_lhs, m_rhs)) {
    std::ostringstream ostr;
    ostr << "The two workspaces are not compatible for algorithm "
         << this->name();
    g_log.error() << ostr.str() << std::endl;
    throw std::invalid_argument(ostr.str());
  }

  // Is the output going to be an EventWorkspace?
  if (m_keepEventWorkspace) {
    // The output WILL be EventWorkspace (this implies lhs is EW or rhs is EW +
    // it gets flipped)
    if (!m_elhs)
      throw std::runtime_error("BinaryOperation:: the output was set to be an "
                               "EventWorkspace (m_keepEventWorkspace == true), "
                               "but the lhs is not an EventWorkspace. There "
                               "must be a mistake in the algorithm. Contact "
                               "the developers.");

    if (m_out == m_lhs) {
      // Will be modifying the EventWorkspace in-place on the lhs. Good.
      if (!m_eout)
        throw std::runtime_error(
            "BinaryOperation:: the output was set to be lhs, and to be an "
            "EventWorkspace (m_keepEventWorkspace == true), but the output is "
            "not an EventWorkspace. There must be a mistake in the algorithm. "
            "Contact the developers.");
    } else {
      // You HAVE to copy the data from lhs to to the output!

      // Create a copy of the lhs workspace
      m_eout = boost::dynamic_pointer_cast<EventWorkspace>(
          API::WorkspaceFactory::Instance().create(
              "EventWorkspace", m_elhs->getNumberHistograms(), 2, 1));
      // Copy geometry, spectra map, etc. over.
      API::WorkspaceFactory::Instance().initializeFromParent(m_elhs, m_eout,
                                                             false);
      // And we copy all the events from the lhs
      m_eout->copyDataFrom(*m_elhs);
      // Make sure m_out still points to the same as m_eout;
      m_out = boost::dynamic_pointer_cast<API::MatrixWorkspace>(m_eout);
    }

    // Always clear the MRUs.
    m_eout->clearMRU();
    if (m_elhs)
      m_elhs->clearMRU();
    if (m_erhs)
      m_erhs->clearMRU();

  } else {
    // ---- Output will be WS2D -------

    // We need to create a new workspace for the output if:
    //   (a) the output workspace hasn't been set to one of the input ones, or
    //   (b) it has been, but it's not the correct dimensions
    if ((m_out != m_lhs && m_out != m_rhs) ||
        (m_out == m_rhs && (m_lhs->size() > m_rhs->size()))) {
      // Make sure to delete anything that might be in the output name.
      // Removed ahead of 2.0 release to avoid problems detailed in trac #4630.
      // Hopefully temporary (see #4635).
      //          if
      //          (AnalysisDataService::Instance().doesExist(getPropertyValue(outputPropName()
      //          )))
      //            AnalysisDataService::Instance().remove(getPropertyValue(outputPropName()
      //            ));
      m_out = WorkspaceFactory::Instance().create(m_lhs);
    }
  }

  // only overridden for some operations (plus and minus at the time of writing)
  operateOnRun(m_lhs->run(), m_rhs->run(), m_out->mutableRun());

  // Initialise the progress reporting object
  m_progress = new Progress(this, 0.0, 1.0, m_lhs->getNumberHistograms());

  // There are now 4 possible scenarios, shown schematically here:
  // xxx x   xxx xxx   xxx xxx   xxx x
  // xxx   , xxx xxx , xxx     , xxx x
  // xxx   , xxx xxx   xxx       xxx x
  // So work out which one we have and call the appropriate function

  // Single value workspace on the right : if it is an EventWorkspace with 1
  // spectrum, 1 bin, it is treated as a scalar
  if ((m_rhs->size() == 1) && !m_do2D_even_for_SingleColumn_on_rhs) {
    doSingleValue();                            // m_lhs,m_rhs,m_out
  } else if (m_rhs->getNumberHistograms() == 1) // Single spectrum on rhs
  {
    doSingleSpectrum();
  }
  // Single column on rhs; if the RHS is an event workspace with one bin, it is
  // treated as a scalar.
  else if ((m_rhs->blocksize() == 1) && !m_do2D_even_for_SingleColumn_on_rhs) {
    m_indicesToMask.reserve(m_out->getNumberHistograms());
    doSingleColumn();
  } else // The two are both 2D and should be the same size (except if LHS is an
         // event workspace)
  {
    m_indicesToMask.reserve(m_out->getNumberHistograms());

    bool mismatchedSpectra =
        (m_AllowDifferentNumberSpectra &&
         (m_rhs->getNumberHistograms() != m_lhs->getNumberHistograms()));
    do2D(mismatchedSpectra);
  }

  applyMaskingToOutput(m_out);
  setOutputUnits(m_lhs, m_rhs, m_out);

  // Assign the result to the output workspace property
  setProperty(outputPropName(), m_out);

  return;
}

//--------------------------------------------------------------------------------------------
/**
 * Execute a binary operation on events. Should be overridden.
 * @param lhs :: left-hand event workspace
 * @param rhs :: right-hand event workspace
 */
void BinaryOperation::execEvent(DataObjects::EventWorkspace_const_sptr lhs,
                                DataObjects::EventWorkspace_const_sptr rhs) {
  UNUSED_ARG(lhs);
  UNUSED_ARG(rhs);
  // This should never happen
  throw Exception::NotImplementedError(
      "BinaryOperation::execEvent() is not implemented for this operation.");
}

//--------------------------------------------------------------------------------------------
/**
 * Return true if the two workspaces are compatible for this operation
 * Virtual: will be overridden as needed.
 * @param lhs :: left-hand workspace to check
 * @param rhs :: right-hand workspace to check
 * @return flag for the compatibility to the two workspaces
 */
bool BinaryOperation::checkCompatibility(
    const API::MatrixWorkspace_const_sptr lhs,
    const API::MatrixWorkspace_const_sptr rhs) const {
  Unit_const_sptr lhs_unit;
  Unit_const_sptr rhs_unit;
  if (lhs->axes() && rhs->axes()) // If one of these is a WorkspaceSingleValue
                                  // then we don't want to check units match
  {
    lhs_unit = lhs->getAxis(0)->unit();
    rhs_unit = rhs->getAxis(0)->unit();
  }

  const std::string lhs_unitID = (lhs_unit ? lhs_unit->unitID() : "");
  const std::string rhs_unitID = (rhs_unit ? rhs_unit->unitID() : "");

  // Check the workspaces have the same units and distribution flag
  if (lhs_unitID != rhs_unitID && lhs->blocksize() > 1 &&
      rhs->blocksize() > 1) {
    g_log.error("The two workspace are not compatible because they have "
                "different units on the X axis.");
    return false;
  }

  // Check the size compatibility
  std::string checkSizeCompatibilityResult = checkSizeCompatibility(lhs, rhs);
  if (!checkSizeCompatibilityResult.empty()) {
    throw std::invalid_argument(checkSizeCompatibilityResult);
  }

  return true;
}

//--------------------------------------------------------------------------------------------
/** Return true if the two workspaces can be treated as event workspaces
 * for the binary operation. If so, execEvent() will be called.
 * (e.g. Plus algorithm will concatenate event lists)
 * @param lhs :: left-hand event workspace to check
 * @param rhs :: right-hand event workspace to check
 * @return false by default; will be overridden by specific algorithms
 */
bool BinaryOperation::checkEventCompatibility(
    const API::MatrixWorkspace_const_sptr lhs,
    const API::MatrixWorkspace_const_sptr rhs) {
  UNUSED_ARG(lhs);
  UNUSED_ARG(rhs);
  return false;
}

//--------------------------------------------------------------------------------------------
/** Performs a simple check to see if the sizes of two workspaces are compatible
 * for a binary operation
 *  In order to be size compatible then the larger workspace
 *  must divide be the size of the smaller workspace leaving no remainder
 *  @param lhs :: the first workspace to compare
 *  @param rhs :: the second workspace to compare
 *  @retval "" The two workspaces are size compatible
 *  @retval "<reason why not compatible>" The two workspaces are NOT size
 * compatible
 */
std::string BinaryOperation::checkSizeCompatibility(
    const API::MatrixWorkspace_const_sptr lhs,
    const API::MatrixWorkspace_const_sptr rhs) const {
  const size_t lhsSize = lhs->size();
  const size_t rhsSize = rhs->size();
  // A SingleValueWorkspace on the right matches anything
  if (rhsSize == 1)
    return "";
  // The lhs must not be smaller than the rhs
  if (lhsSize < rhsSize)
    return "Left hand side smaller than right hand side.";

  // Did checkRequirements() tell us that the X histogram size did not matter?
  if (!m_matchXSize) {
    // If so, only the vertical # needs to match

    if (lhs->getNumberHistograms() == rhs->getNumberHistograms()) {
      return "";
    } else {
      return "Number of histograms not identical.";
    }
  }
  // Otherwise they must match both ways, or horizontally or vertically with the
  // other rhs dimension=1
  if (rhs->blocksize() == 1 &&
      lhs->getNumberHistograms() == rhs->getNumberHistograms())
    return "";
  // Past this point, we require the X arrays to match. Note this only checks
  // the first spectrum
  if (!WorkspaceHelpers::matchingBins(lhs, rhs, true)) {
    return "X arrays must match when performing this operation on a 2D "
           "workspaces.";
  }

  const size_t rhsSpec = rhs->getNumberHistograms();

  if (lhs->blocksize() == rhs->blocksize()) {
    if (rhsSpec == 1 || lhs->getNumberHistograms() == rhsSpec) {
      return "";
    } else {
      // can't be more specific as if this is reached both failed and only one
      // or both are needed
      return "Left and right sides should contain the same amount of spectra "
             "or the right side should contian only one spectra.";
    }
  } else {
    // blocksize check failed, but still check the number of spectra to see if
    // that was wrong too
    if (rhsSpec == 1 || lhs->getNumberHistograms() == rhsSpec) {
      return "Number of y values not equal on left and right sides.";
    } else {
      // can't be more specific as if this is reached both failed and only one
      // or both are needed
      return "Number of y values not equal on left and right sides and the "
             "right side contained neither only one spectra or the same amount "
             "of spectra as the left.";
    }
  }
}

//--------------------------------------------------------------------------------------------
/**
 * Checks if the spectra at the given index of either input workspace is masked.
 * If so then the output spectra has zeroed data
 * and is also masked.
 * @param lhs :: A pointer to the left-hand operand
 * @param rhs :: A pointer to the right-hand operand
 * @param index :: The workspace index to check
 * @param out :: A pointer to the output workspace
 * @returns True if further processing is not required on the spectra, false if
 * the binary operation should be performed.
 */
bool
BinaryOperation::propagateSpectraMask(const API::MatrixWorkspace_const_sptr lhs,
                                      const API::MatrixWorkspace_const_sptr rhs,
                                      const int64_t index,
                                      API::MatrixWorkspace_sptr out) {
  bool continueOp(true);
  IDetector_const_sptr det_lhs, det_rhs;
  try {
    det_lhs = lhs->getDetector(index);
    det_rhs = rhs->getDetector(index);
  } catch (std::runtime_error &) {
  } catch (std::domain_error &) {
    // try statement will throw a domain_error when the axis is not a spectra
    // axis.
    return continueOp;
  }
  if ((det_lhs && det_lhs->isMasked()) || (det_rhs && det_rhs->isMasked())) {
    continueOp = false;
    out->maskWorkspaceIndex(index);
  }
  return continueOp;
}

/**
 * Called when the rhs operand is a single value.
 *  Loops over the lhs workspace calling the abstract binary operation function
 * with a single number as the rhs operand.
 */
void BinaryOperation::doSingleValue() {
  // Don't propate masking from the rhs here - it would be decidedly odd if the
  // single value was masked

  // Pull out the single value and its error
  const double rhsY = m_rhs->readY(0)[0];
  const double rhsE = m_rhs->readE(0)[0];

  // Now loop over the spectra of the left hand side calling the virtual
  // function
  const int64_t numHists = m_lhs->getNumberHistograms();

  if (m_eout) {
    // ---- The output is an EventWorkspace ------
    PARALLEL_FOR3(m_lhs, m_rhs, m_out)
    for (int64_t i = 0; i < numHists; ++i) {
      PARALLEL_START_INTERUPT_REGION
      m_out->setX(i, m_lhs->refX(i));
      performEventBinaryOperation(m_eout->getEventList(i), rhsY, rhsE);
      m_progress->report(this->name());
      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION
  } else {
    // ---- Histogram Output -----
    PARALLEL_FOR3(m_lhs, m_rhs, m_out)
    for (int64_t i = 0; i < numHists; ++i) {
      PARALLEL_START_INTERUPT_REGION
      m_out->setX(i, m_lhs->refX(i));
      // Get reference to output vectors here to break any sharing outside the
      // function call below
      // where the order of argument evaluation is not guaranteed (if it's L->R
      // there would be a data race)
      MantidVec &outY = m_out->dataY(i);
      MantidVec &outE = m_out->dataE(i);
      performBinaryOperation(m_lhs->readX(i), m_lhs->readY(i), m_lhs->readE(i),
                             rhsY, rhsE, outY, outE);
      m_progress->report(this->name());
      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION
  }
}

//--------------------------------------------------------------------------------------------
/** Called when the m_rhs operand is a 2D workspace of single values.
 *  Loops over the workspaces calling the abstract binary operation function
 * with a single number as the m_rhs operand.
 */
void BinaryOperation::doSingleColumn() {
  // Don't propate masking from the m_rhs here - it would be decidedly odd if
  // the single bin was masked

  // Now loop over the spectra of the left hand side pulling m_out the single
  // value from each m_rhs 'spectrum'
  // and then calling the virtual function
  const int64_t numHists = m_lhs->getNumberHistograms();
  if (m_eout) {
    // ---- The output is an EventWorkspace ------
    PARALLEL_FOR3(m_lhs, m_rhs, m_out)
    for (int64_t i = 0; i < numHists; ++i) {
      PARALLEL_START_INTERUPT_REGION
      const double rhsY = m_rhs->readY(i)[0];
      const double rhsE = m_rhs->readE(i)[0];

      // m_out->setX(i, m_lhs->refX(i)); //unnecessary - that was copied before.
      if (propagateSpectraMask(m_lhs, m_rhs, i, m_out)) {
        performEventBinaryOperation(m_eout->getEventList(i), rhsY, rhsE);
      }
      m_progress->report(this->name());
      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION
  } else {
    // ---- Histogram Output -----
    PARALLEL_FOR3(m_lhs, m_rhs, m_out)
    for (int64_t i = 0; i < numHists; ++i) {
      PARALLEL_START_INTERUPT_REGION
      const double rhsY = m_rhs->readY(i)[0];
      const double rhsE = m_rhs->readE(i)[0];

      m_out->setX(i, m_lhs->refX(i));
      if (propagateSpectraMask(m_lhs, m_rhs, i, m_out)) {
        // Get reference to output vectors here to break any sharing outside the
        // function call below
        // where the order of argument evaluation is not guaranteed (if it's
        // L->R there would be a data race)
        MantidVec &outY = m_out->dataY(i);
        MantidVec &outE = m_out->dataE(i);
        performBinaryOperation(m_lhs->readX(i), m_lhs->readY(i),
                               m_lhs->readE(i), rhsY, rhsE, outY, outE);
      }
      m_progress->report(this->name());
      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION
  }
}

//--------------------------------------------------------------------------------------------
/** Called when the m_rhs operand is a single spectrum.
 *  Loops over the lhs workspace calling the abstract binary operation function.
 */
void BinaryOperation::doSingleSpectrum() {

  // Propagate any masking first or it could mess up the numbers
  // TODO: Check if this works for event workspaces...
  propagateBinMasks(m_rhs, m_out);

  if (m_eout) {
    // ----------- The output is an EventWorkspace -------------

    if (m_erhs) {
      // -------- The rhs is ALSO an EventWorkspace --------

      // Pull out the single eventList on the right
      const EventList &rhs_spectrum = m_erhs->getEventList(0);

      // Now loop over the spectra of the left hand side calling the virtual
      // function
      const int64_t numHists = m_lhs->getNumberHistograms();
      PARALLEL_FOR3(m_lhs, m_rhs, m_out)
      for (int64_t i = 0; i < numHists; ++i) {
        PARALLEL_START_INTERUPT_REGION
        // m_out->setX(i,m_lhs->refX(i)); //unnecessary - that was copied
        // before.

        // Perform the operation on the event list on the output (== lhs)
        performEventBinaryOperation(m_eout->getEventList(i), rhs_spectrum);
        m_progress->report(this->name());
        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION
    } else {
      // -------- The rhs is a histogram ---------
      // Pull m_out the m_rhs spectrum
      const MantidVec &rhsX = m_rhs->readX(0);
      const MantidVec &rhsY = m_rhs->readY(0);
      const MantidVec &rhsE = m_rhs->readE(0);

      // Now loop over the spectra of the left hand side calling the virtual
      // function
      const int64_t numHists = m_lhs->getNumberHistograms();

      PARALLEL_FOR3(m_lhs, m_rhs, m_out)
      for (int64_t i = 0; i < numHists; ++i) {
        PARALLEL_START_INTERUPT_REGION
        // m_out->setX(i,m_lhs->refX(i)); //unnecessary - that was copied
        // before.
        // Perform the operation on the event list on the output (== lhs)
        performEventBinaryOperation(m_eout->getEventList(i), rhsX, rhsY, rhsE);
        m_progress->report(this->name());
        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION
    }

  } else {
    // -------- The output is a histogram ----------
    // (inputs can be EventWorkspaces, but their histogram representation
    //  will be used instead)

    // Pull m_out the m_rhs spectrum
    const MantidVec &rhsY = m_rhs->readY(0);
    const MantidVec &rhsE = m_rhs->readE(0);

    // Now loop over the spectra of the left hand side calling the virtual
    // function
    const int64_t numHists = m_lhs->getNumberHistograms();

    PARALLEL_FOR3(m_lhs, m_rhs, m_out)
    for (int64_t i = 0; i < numHists; ++i) {
      PARALLEL_START_INTERUPT_REGION
      m_out->setX(i, m_lhs->refX(i));
      // Get reference to output vectors here to break any sharing outside the
      // function call below
      // where the order of argument evaluation is not guaranteed (if it's L->R
      // there would be a data race)
      MantidVec &outY = m_out->dataY(i);
      MantidVec &outE = m_out->dataE(i);
      performBinaryOperation(m_lhs->readX(i), m_lhs->readY(i), m_lhs->readE(i),
                             rhsY, rhsE, outY, outE);
      m_progress->report(this->name());
      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION
  }
}

//--------------------------------------------------------------------------------------------
/** Called when the two workspaces are the same size.
 *  Loops over the workspaces extracting the appropriate spectra and calling the
 *abstract binary operation function.
 *
 *  @param mismatchedSpectra :: allow the # of spectra to not be the same. Will
 *use the
 *      detector IDs to find the corresponding spectrum on RHS
 */
void BinaryOperation::do2D(bool mismatchedSpectra) {
  BinaryOperationTable_sptr table;
  if (mismatchedSpectra) {
    table = BinaryOperation::buildBinaryOperationTable(m_lhs, m_rhs);
  }

  // Propagate any masking first or it could mess up the numbers
  // TODO: Check if this works for event workspaces...
  propagateBinMasks(m_rhs, m_out);

  if (m_eout) {
    // ----------- The output is an EventWorkspace -------------

    if (m_erhs && !m_useHistogramForRhsEventWorkspace) {
      // ------------ The rhs is ALSO an EventWorkspace ---------------
      // Now loop over the spectra of each one calling the virtual function
      const int64_t numHists = m_lhs->getNumberHistograms();
      PARALLEL_FOR3(m_lhs, m_rhs, m_out)
      for (int64_t i = 0; i < numHists; ++i) {
        PARALLEL_START_INTERUPT_REGION
        m_progress->report(this->name());

        int64_t rhs_wi = i;
        if (mismatchedSpectra && table) {
          rhs_wi = (*table)[i];
          if (rhs_wi < 0)
            continue;
        } else {
          // Check for masking except when mismatched sizes
          if (!propagateSpectraMask(m_lhs, m_rhs, i, m_out))
            continue;
        }
        // Reach here? Do the division
        // Perform the operation on the event list on the output (== lhs)
        performEventBinaryOperation(m_eout->getEventList(i),
                                    m_erhs->getEventList(rhs_wi));

        // Free up memory on the RHS if that is possible
        if (m_ClearRHSWorkspace)
          const_cast<EventList &>(m_erhs->getEventList(rhs_wi)).clear();
        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION
    } else {
      // -------- The rhs is a histogram, or we want to use the histogram
      // representation of it ---------

      // Now loop over the spectra of each one calling the virtual function
      const int64_t numHists = m_lhs->getNumberHistograms();

      PARALLEL_FOR3(m_lhs, m_rhs, m_out)
      for (int64_t i = 0; i < numHists; ++i) {
        PARALLEL_START_INTERUPT_REGION
        m_progress->report(this->name());
        int64_t rhs_wi = i;
        if (mismatchedSpectra && table) {
          rhs_wi = (*table)[i];
          if (rhs_wi < 0)
            continue;
        } else {
          // Check for masking except when mismatched sizes
          if (!propagateSpectraMask(m_lhs, m_rhs, i, m_out))
            continue;
        }

        // Reach here? Do the division
        performEventBinaryOperation(m_eout->getEventList(i),
                                    m_rhs->readX(rhs_wi), m_rhs->readY(rhs_wi),
                                    m_rhs->readE(rhs_wi));

        // Free up memory on the RHS if that is possible
        if (m_ClearRHSWorkspace)
          const_cast<EventList &>(m_erhs->getEventList(rhs_wi)).clear();

        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION
    }

  } else {
    // -------- The output is a histogram ----------
    // (inputs can be EventWorkspaces, but their histogram representation
    //  will be used instead)

    // Now loop over the spectra of each one calling the virtual function
    const int64_t numHists = m_lhs->getNumberHistograms();

    PARALLEL_FOR3(m_lhs, m_rhs, m_out)
    for (int64_t i = 0; i < numHists; ++i) {
      PARALLEL_START_INTERUPT_REGION
      m_progress->report(this->name());
      m_out->setX(i, m_lhs->refX(i));
      int64_t rhs_wi = i;
      if (mismatchedSpectra && table) {
        rhs_wi = (*table)[i];
        if (rhs_wi < 0)
          continue;
      } else {
        // Check for masking except when mismatched sizes
        if (!propagateSpectraMask(m_lhs, m_rhs, i, m_out))
          continue;
      }
      // Reach here? Do the division
      // Get reference to output vectors here to break any sharing outside the
      // function call below
      // where the order of argument evaluation is not guaranteed (if it's L->R
      // there would be a data race)
      MantidVec &outY = m_out->dataY(i);
      MantidVec &outE = m_out->dataE(i);
      performBinaryOperation(m_lhs->readX(i), m_lhs->readY(i), m_lhs->readE(i),
                             m_rhs->readY(rhs_wi), m_rhs->readE(rhs_wi), outY,
                             outE);

      // Free up memory on the RHS if that is possible
      if (m_ClearRHSWorkspace)
        const_cast<EventList &>(m_erhs->getEventList(rhs_wi)).clear();

      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION
  }
  // Make sure we don't use outdated MRU
  if (m_ClearRHSWorkspace)
    m_erhs->clearMRU();
}

/** Copies any bin masking from the smaller/rhs input workspace to the output.
 *  Masks on the other input workspace are copied automatically by the workspace
 * factory.
 *  @param rhs :: The workspace which is the right hand operand
 *  @param out :: The result workspace
 */
void
BinaryOperation::propagateBinMasks(const API::MatrixWorkspace_const_sptr rhs,
                                   API::MatrixWorkspace_sptr out) {
  const int64_t outHists = out->getNumberHistograms();
  const int64_t rhsHists = rhs->getNumberHistograms();
  for (int64_t i = 0; i < outHists; ++i) {
    // Copy over masks from the rhs, if any exist.
    // If rhs is single spectrum, copy masks from that to all spectra in the
    // output.
    if (rhs->hasMaskedBins((rhsHists == 1) ? 0 : i)) {
      const MatrixWorkspace::MaskList &masks =
          rhs->maskedBins((rhsHists == 1) ? 0 : i);
      MatrixWorkspace::MaskList::const_iterator it;
      for (it = masks.begin(); it != masks.end(); ++it) {
        out->flagMasked(i, it->first, it->second);
      }
    }
  }
}

//---------------------------------------------------------------------------------------------
/**
 * Apply the requested masking to the output workspace
 * @param out :: The workspace to mask
 */
void BinaryOperation::applyMaskingToOutput(API::MatrixWorkspace_sptr out) {
  int64_t nindices = static_cast<int64_t>(m_indicesToMask.size());
  ParameterMap &pmap = out->instrumentParameters();
  PARALLEL_FOR1(out)
  for (int64_t i = 0; i < nindices; ++i) {
    if (!m_parallelException && !m_cancel) {
      try {
        IDetector_const_sptr det_out = out->getDetector(m_indicesToMask[i]);
        PARALLEL_CRITICAL(BinaryOperation_masking) {
          pmap.addBool(det_out.get(), "masked", true);
        }
      } /* End of try block in PARALLEL_START_INTERUPT_REGION */
      catch (Kernel::Exception::NotFoundError) { // detector not found, do
                                                 // nothing, go further
      } catch (std::runtime_error &ex) {
        if (!m_parallelException) {
          m_parallelException = true;
          g_log.error() << this->name() << ": " << ex.what() << "\n";
        }
      } catch (...) {
        m_parallelException = true;
      }

    } // End of if block in PARALLEL_START_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION
}

// ------- Default implementations of Event binary operations --------

/**
 * Carries out the binary operation IN-PLACE on a single EventList,
 * with another EventList as the right-hand operand.
 * The event lists simply get appended.
 *
 *  @param lhs :: Reference to the EventList that will be modified in place.
 *  @param rhs :: Const reference to the EventList on the right hand side.
 */
void BinaryOperation::performEventBinaryOperation(
    DataObjects::EventList &lhs, const DataObjects::EventList &rhs) {
  UNUSED_ARG(lhs);
  UNUSED_ARG(rhs);
  throw Exception::NotImplementedError(
      "BinaryOperation::performEventBinaryOperation() not implemented.");
}

/**
 * Carries out the binary operation IN-PLACE on a single EventList,
 * with another (histogrammed) spectrum as the right-hand operand.
 *
 *  @param lhs :: Reference to the EventList that will be modified in place.
 *  @param rhsX :: The vector of rhs X bin boundaries
 *  @param rhsY :: The vector of rhs data values
 *  @param rhsE :: The vector of rhs error values
 */
void BinaryOperation::performEventBinaryOperation(DataObjects::EventList &lhs,
                                                  const MantidVec &rhsX,
                                                  const MantidVec &rhsY,
                                                  const MantidVec &rhsE) {
  UNUSED_ARG(lhs);
  UNUSED_ARG(rhsX);
  UNUSED_ARG(rhsY);
  UNUSED_ARG(rhsE);
  throw Exception::NotImplementedError(
      "BinaryOperation::performEventBinaryOperation() not implemented.");
}

/**
 * Carries out the binary operation IN-PLACE on a single EventList,
 * with a single (double) value as the right-hand operand
 *
 *  @param lhs :: Reference to the EventList that will be modified in place.
 *  @param rhsY :: The rhs data value
 *  @param rhsE :: The rhs error value
 */
void BinaryOperation::performEventBinaryOperation(DataObjects::EventList &lhs,
                                                  const double &rhsY,
                                                  const double &rhsE) {
  UNUSED_ARG(lhs);
  UNUSED_ARG(rhsY);
  UNUSED_ARG(rhsE);
  throw Exception::NotImplementedError(
      "BinaryOperation::performEventBinaryOperation() not implemented.");
}

//---------------------------------------------------------------------------------------------
/**
 * Get the type of operand from a workspace
 * @param ws :: workspace to check
 * @return OperandType describing what type of workspace it will be operated as.
 */
OperandType
BinaryOperation::getOperandType(const API::MatrixWorkspace_const_sptr ws) {
  // An event workspace?
  EventWorkspace_const_sptr ews =
      boost::dynamic_pointer_cast<const EventWorkspace>(ws);
  if (ews)
    return eEventList;

  // If the workspace has no axes, then it is a WorkspaceSingleValue
  if (!ws->axes())
    return eNumber;

  // TODO: Check if it is a single-colum one, then
  //  return Number;

  // Otherwise, we take it as a histogram (workspace 2D)
  return eHistogram;
}

//---------------------------------------------------------------------------------------------
/** Check what operation will be needed in order to apply the operation
 * to these two types of workspaces. This function must be overridden
 * and checked against all 9 possible combinations.
 *
 * Must set: m_matchXSize, m_flipSides, m_keepEventWorkspace
 */
void BinaryOperation::checkRequirements() {

  // In general, the X sizes have to match.
  //  (only for some EventWorkspace cases does it NOT matter...)
  m_matchXSize = true;

  // Operations are not always commutative! Don't flip sides.
  m_flipSides = false;

  // And in general, EventWorkspaces get turned to Workspace2D
  m_keepEventWorkspace = false;

  // This will be set to true for Divide/Multiply
  m_useHistogramForRhsEventWorkspace = false;
}

//---------------------------------------------------------------------------------------------
/** Build up an BinaryOperationTable for performing a binary operation
 * e.g. lhs = (lhs + rhs)
 * where the spectra in rhs are to go into lhs.
 * This function looks to match the detector IDs in rhs to those in the lhs.
 *
 * @param lhs :: matrix workspace in which the operation is being done.
 * @param rhs :: matrix workspace on the right hand side of the operand
 * @return map from detector ID to workspace index for the RHS workspace.
 *        NULL if there is not a 1:1 mapping from detector ID to workspace index
 *(e.g more than one detector per pixel).
 */
BinaryOperation::BinaryOperationTable_sptr
BinaryOperation::buildBinaryOperationTable(
    const MatrixWorkspace_const_sptr &lhs,
    const MatrixWorkspace_const_sptr &rhs) {
  // An addition table is a list of pairs:
  //  First int = workspace index in the EW being added
  //  Second int = workspace index to which it will be added in the OUTPUT EW.
  //  -1 if it should add a new entry at the end.
  auto table = boost::make_shared<BinaryOperationTable>();

  int rhs_nhist = static_cast<int>(rhs->getNumberHistograms());
  int lhs_nhist = static_cast<int>(lhs->getNumberHistograms());

  // Initialize the table; filled with -1 meaning no match
  table->resize(lhs_nhist, -1);

  const detid2index_map rhs_det_to_wi = rhs->getDetectorIDToWorkspaceIndexMap();

  PARALLEL_FOR_NO_WSP_CHECK()
  for (int lhsWI = 0; lhsWI < lhs_nhist; lhsWI++) {
    bool done = false;

    // List of detectors on lhs side
    const std::set<detid_t> &lhsDets =
        lhs->getSpectrum(lhsWI)->getDetectorIDs();

    // ----------------- Matching Workspace Indices and Detector IDs
    // --------------------------------------
    // First off, try to match the workspace indices. Most times, this will be
    // ok right away.
    int64_t rhsWI = lhsWI;
    if (rhsWI < rhs_nhist) // don't go out of bounds
    {
      // Get the detector IDs at that workspace index.
      const std::set<detid_t> &rhsDets =
          rhs->getSpectrum(rhsWI)->getDetectorIDs();

      // Checks that lhsDets is a subset of rhsDets
      if (std::includes(rhsDets.begin(), rhsDets.end(), lhsDets.begin(),
                        lhsDets.end())) {
        // We found the workspace index right away. No need to keep looking
        (*table)[lhsWI] = rhsWI;
        done = true;
      }
    }

    // ----------------- Scrambled Detector IDs with one Detector per Spectrum
    // --------------------------------------
    if (!done && (lhsDets.size() == 1)) {
      // Didn't find it. Try to use the RHS map.

      // First, we have to get the (single) detector ID of the LHS
      std::set<detid_t>::const_iterator lhsDets_it = lhsDets.begin();
      detid_t lhs_detector_ID = *lhsDets_it;

      // Now we use the RHS map to find it. This only works if both the lhs and
      // rhs have 1 detector per pixel
      detid2index_map::const_iterator map_it =
          rhs_det_to_wi.find(lhs_detector_ID);
      if (map_it != rhs_det_to_wi.end()) {
        rhsWI = map_it->second; // This is the workspace index in the RHS that
                                // matched lhs_detector_ID
      } else {
        // Did not find it!
        rhsWI = -1; // Marker to mean its not in the LHS.

        //            std::ostringstream mess;
        //            mess << "BinaryOperation: cannot find a RHS spectrum that
        //            contains the detectors in LHS workspace index " << lhsWI
        //            << "\n";
        //            throw std::runtime_error(mess.str());
      }
      (*table)[lhsWI] = rhsWI;
      done = true; // Great, we did it.
    }

    // ----------------- LHS detectors are subset of RHS, which are Grouped
    // --------------------------------------
    if (!done) {

      // Didn't find it? Now we need to iterate through the output workspace to
      //  match the detector ID.
      // NOTE: This can be SUPER SLOW!
      for (rhsWI = 0; rhsWI < static_cast<int64_t>(rhs_nhist); rhsWI++) {
        const std::set<detid_t> &rhsDets =
            rhs->getSpectrum(rhsWI)->getDetectorIDs();

        // Checks that lhsDets is a subset of rhsDets
        if (std::includes(rhsDets.begin(), rhsDets.end(), lhsDets.begin(),
                          lhsDets.end())) {
          // This one is right. Now we can stop looking.
          (*table)[lhsWI] = rhsWI;
          done = true;
          continue;
        }
      }
    }

    // ------- Still nothing ! -----------
    if (!done) {
      (*table)[lhsWI] = -1;

      //          std::ostringstream mess;
      //          mess << "BinaryOperation: cannot find a RHS spectrum that
      //          contains the detectors in LHS workspace index " << lhsWI <<
      //          "\n";
      //          throw std::runtime_error(mess.str());
    }
  }

  return table;
}

} // namespace Algorithms
} // namespace Mantid