ConvertAxisByFormula.cpp

#include "MantidAlgorithms/ConvertAxisByFormula.h"
#include "MantidAPI/CommonBinsValidator.h"
#include "MantidAPI/RefAxis.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidGeometry/muParser_Silent.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/UnitFactory.h"

#include <boost/shared_ptr.hpp>
#include <boost/make_shared.hpp>
#include <sstream>
#include <algorithm>

namespace Mantid {
namespace Algorithms {

using namespace Kernel;
using namespace API;

// Register the class into the algorithm factory
DECLARE_ALGORITHM(ConvertAxisByFormula)

const std::string ConvertAxisByFormula::name() const {
  return ("ConvertAxisByFormula");
}

int ConvertAxisByFormula::version() const { return (1); }

const std::string ConvertAxisByFormula::category() const {
  return "Transforms\\Axes";
}

/** Initialisation method. Declares properties to be used in algorithm.
*
*/
void ConvertAxisByFormula::init() {
  declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
                      "InputWorkspace", "", Direction::Input),
                  "Name of the input workspace");
  declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
                      "OutputWorkspace", "", Direction::Output),
                  "Name of the output workspace");

  std::vector<std::string> axisOptions;
  axisOptions.emplace_back("X");
  axisOptions.emplace_back("Y");
  declareProperty("Axis", "X",
                  boost::make_shared<StringListValidator>(axisOptions),
                  "The axis to modify");

  declareProperty("Formula", "", "The formula to use to convert the values, x "
                                 "or y may be used to refer to the axis "
                                 "values.  l1, l2, twotheta and signedtwotheta"
                                 "may be used to provide values from the "
                                 "instrument geometry.");
  declareProperty(
      "AxisTitle", "",
      "The label of he new axis. If not set then the title will not change.");
  declareProperty(
      "AxisUnits", "",
      "The units of the new axis. If not set then the unit will not change");
}

/** Execution of the algorithm
*
*/
void ConvertAxisByFormula::exec() {
  // get the property values
  MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
  std::string axis = getProperty("Axis");
  std::string formula = getProperty("Formula");
  std::string axisTitle = getProperty("AxisTitle");
  std::string axisUnits = getProperty("AxisUnits");

  // Just overwrite if the change is in place
  MatrixWorkspace_sptr outputWs = getProperty("OutputWorkspace");
  if (outputWs != inputWs) {
    IAlgorithm_sptr duplicate =
        createChildAlgorithm("CloneWorkspace", 0.0, 0.6);
    duplicate->initialize();
    duplicate->setProperty<Workspace_sptr>(
        "InputWorkspace", boost::dynamic_pointer_cast<Workspace>(inputWs));
    duplicate->execute();
    Workspace_sptr temp = duplicate->getProperty("OutputWorkspace");
    outputWs = boost::dynamic_pointer_cast<MatrixWorkspace>(temp);

    setProperty("OutputWorkspace", outputWs);
  }

  // Get the axes - assume X
  Axis *axisPtr = outputWs->getAxis(0);
  Axis *otherAxisPtr = outputWs->getAxis(1);
  if (axis == "Y") {
    std::swap(axisPtr, otherAxisPtr);
  }

  if (!axisPtr->isNumeric()) {
    throw std::invalid_argument("This algorithm only operates on numeric axes");
  }

  bool isRaggedBins = false;
  bool isRefAxis = false;
  RefAxis *refAxisPtr = dynamic_cast<RefAxis *>(axisPtr);
  if (refAxisPtr != nullptr) {
    CommonBinsValidator sameBins;
    if (!sameBins.isValid(outputWs).empty()) {
      isRaggedBins = true;
    }
    isRefAxis = true;
  }

  // validate - if formula contains geometry variables then the other axis must
  // be a spectraAxis
  SpectraAxis *spectrumAxisPtr = dynamic_cast<SpectraAxis *>(otherAxisPtr);
  std::vector<Variable_ptr> variables;
  variables.reserve(8);
  // axis value lookups
  variables.push_back(boost::make_shared<Variable>("x", false));
  variables.push_back(boost::make_shared<Variable>("X", false));
  variables.push_back(boost::make_shared<Variable>("y", false));
  variables.push_back(boost::make_shared<Variable>("Y", false));
  // geometry lookups
  variables.push_back(boost::make_shared<Variable>("twotheta", true));
  variables.push_back(boost::make_shared<Variable>("signedtwotheta", true));
  variables.push_back(boost::make_shared<Variable>("l1", true));
  variables.push_back(boost::make_shared<Variable>("l2", true));

  bool isGeometryRequired = false;
  for (auto variablesIter = variables.begin();
       variablesIter != variables.end();) {
    if (formula.find((*variablesIter)->name) != std::string::npos) {
      if ((*variablesIter)->isGeometric) {
        if (!spectrumAxisPtr) {
          throw std::invalid_argument("To use geometry values in the equation, "
                                      "the axis not being converted must be a "
                                      "Spectra Axis.");
        }
        isGeometryRequired = true;
      }
      ++variablesIter;
    } else {
      // remove those that are not used
      variablesIter = variables.erase(variablesIter);
    }
  }

  // Create muparser
  mu::Parser p;
  try {
    // set parameter lookups for the axis value
    for (const auto &variable : variables) {
      p.DefineVar(variable->name, &(variable->value));
    }
    // set some constants
    double pi = M_PI;
    p.DefineVar("pi", &pi);
    double h = PhysicalConstants::h;
    p.DefineVar("h", &h);
    double hBar = PhysicalConstants::h_bar;
    p.DefineVar("h_bar", &hBar);
    double g = PhysicalConstants::g;
    p.DefineVar("g", &g);
    double mN = PhysicalConstants::NeutronMass;
    p.DefineVar("mN", &mN);
    double mNAMU = PhysicalConstants::NeutronMassAMU;
    p.DefineVar("mNAMU", &mNAMU);

    p.SetExpr(formula);
  } catch (mu::Parser::exception_type &e) {
    std::stringstream ss;
    ss << "Cannot process the formula"
       << ". Muparser error message is: " << e.GetMsg();
    throw std::invalid_argument(ss.str());
  }
  if (isRefAxis) {
    if ((isRaggedBins) || (isGeometryRequired)) {
      // ragged bins or geometry used - we have to calculate for every spectra
      size_t numberOfSpectra_i = outputWs->getNumberHistograms();
      auto &spectrumInfo = outputWs->mutableSpectrumInfo();

      size_t failedDetectorCount = 0;
      Progress prog(this, 0.6, 1.0, numberOfSpectra_i);
      for (size_t i = 0; i < numberOfSpectra_i; ++i) {
        try {
          MantidVec &vec = outputWs->dataX(i);
          setGeometryValues(spectrumInfo, i, variables);
          calculateValues(p, vec, variables);
        } catch (std::runtime_error &)
        // two possible exceptions runtime error and NotFoundError
        // both handled the same way
        {
          // could not find the geometry info for this spectra
          outputWs->getSpectrum(i).clearData();
          spectrumInfo.setMasked(i, true);
          failedDetectorCount++;
        }
        prog.report();
      }
      if (failedDetectorCount != 0) {
        g_log.information()
            << "Unable to calculate sample-detector distance for "
            << failedDetectorCount << " spectra. Masking spectrum.\n";
      }
    } else {
      // common bins - we only have to calculate once

      // Calculate the new (common) X values
      MantidVec &vec = outputWs->dataX(0);
      calculateValues(p, vec, variables);

      // copy xVals to every spectra
      int64_t numberOfSpectra_i = static_cast<int64_t>(
          outputWs->getNumberHistograms()); // cast to make openmp happy
      auto xVals = outputWs->refX(0);
      Progress prog(this, 0.6, 1.0, numberOfSpectra_i);
      PARALLEL_FOR_IF(Kernel::threadSafe(*outputWs))
      for (int64_t j = 1; j < numberOfSpectra_i; ++j) {
        PARALLEL_START_INTERUPT_REGION
        outputWs->setX(j, xVals);
        prog.report();
        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION
    }
  } else {
    size_t axisLength = axisPtr->length();
    for (size_t i = 0; i < axisLength; ++i) {
      setAxisValue(axisPtr->getValue(i), variables);
      axisPtr->setValue(i, evaluateResult(p));
    }
  }

  // If the units conversion has flipped the ascending direction of X, reverse
  // all the vectors
  size_t midSpectra = outputWs->getNumberHistograms() / 2;
  if (!outputWs->dataX(0).empty() &&
      (outputWs->dataX(0).front() > outputWs->dataX(0).back() ||
       outputWs->dataX(midSpectra).front() >
           outputWs->dataX(midSpectra).back())) {
    g_log.information("Reversing data within the workspace to keep the axes in "
                      "increasing order.");
    this->reverse(outputWs);
  }

  // Set the Unit of the Axis
  if ((!axisUnits.empty()) || (!axisTitle.empty())) {
    try {
      axisPtr->unit() = UnitFactory::Instance().create(axisUnits);
    } catch (Exception::NotFoundError &) {
      if (axisTitle.empty()) {
        axisTitle = axisPtr->unit()->caption();
      }
      if (axisUnits.empty()) {
        axisUnits = axisPtr->unit()->label();
      }
      axisPtr->unit() = boost::make_shared<Units::Label>(axisTitle, axisUnits);
    }
  }
}

void ConvertAxisByFormula::setAxisValue(const double &value,
                                        std::vector<Variable_ptr> &variables) {
  for (const auto &variable : variables) {
    if (!variable->isGeometric) {
      variable->value = value;
    }
  }
}

void ConvertAxisByFormula::calculateValues(
    mu::Parser &p, MantidVec &vec, std::vector<Variable_ptr> variables) {
  MantidVec::iterator iter;
  for (iter = vec.begin(); iter != vec.end(); ++iter) {
    setAxisValue(*iter, variables);
    *iter = evaluateResult(p);
  }
}

void ConvertAxisByFormula::setGeometryValues(
    const API::SpectrumInfo &specInfo, const size_t index,
    std::vector<Variable_ptr> &variables) {
  for (const auto &variable : variables) {
    if (variable->isGeometric) {
      if (variable->name == "twotheta") {
        variable->value = specInfo.twoTheta(index);
      } else if (variable->name == "signedtwotheta") {
        variable->value = specInfo.signedTwoTheta(index);
      } else if (variable->name == "l1") {
        variable->value = specInfo.l1();
      } else if (variable->name == "l2") {
        variable->value = specInfo.l2(index);
      }
    }
  }
}

double ConvertAxisByFormula::evaluateResult(mu::Parser &p) {
  double result;
  try {
    result = p.Eval();
  } catch (mu::Parser::exception_type &e) {
    std::stringstream ss;
    ss << "Failed while processing axis values"
       << ". Muparser error message is: " << e.GetMsg();
    throw std::invalid_argument(ss.str());
  }
  return result;
}

} // namespace Algorithms
} // namespace Mantid