BinaryOperation.cpp

#include "MantidAlgorithms/BinaryOperation.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/FrameworkManager.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceOpOverloads.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidDataObjects/EventList.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/WorkspaceSingleValue.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidKernel/Timer.h"
#include "MantidKernel/Unit.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 {
/** Initialisation method.
 *  Defines input and output workspaces
 *
 */
void BinaryOperation::init() {
  declareProperty(
      Kernel::make_unique<WorkspaceProperty<MatrixWorkspace>>(
          inputPropName1(), "", Direction::Input),
      "The name of the input workspace on the left hand side of the operation");
  declareProperty(Kernel::make_unique<WorkspaceProperty<MatrixWorkspace>>(
                      inputPropName2(), "", Direction::Input),
                  "The name of the input workspace on the right hand side of "
                  "the operation");
  declareProperty(Kernel::make_unique<WorkspaceProperty<MatrixWorkspace>>(
                      outputPropName(), "", Direction::Output),
                  "The name to call the output workspace");
  declareProperty(
      make_unique<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(
      make_unique<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() << '\n';
    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.");