#include "MantidAlgorithms/ConvertAxisByFormula.h" #include "MantidAPI/CommonBinsValidator.h" #include "MantidAPI/IncreasingAxisValidator.h" #include "MantidAPI/RefAxis.h" #include "MantidGeometry/muParser_Silent.h" #include "MantidKernel/ListValidator.h" #include "MantidKernel/UnitFactory.h" #include <boost/shared_ptr.hpp> #include <boost/make_shared.hpp> #include <sstream> 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"); 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 axis int axisIndex = 0; // assume X if (axis == "Y") { axisIndex = 1; } Axis *axisPtr = outputWs->getAxis(axisIndex); if (!axisPtr->isNumeric()) { throw std::invalid_argument("This algorithm only operates on numeric axes"); } bool isRefAxis = false; RefAxis *refAxisPtr = dynamic_cast<RefAxis *>(axisPtr); if (refAxisPtr != nullptr) { CommonBinsValidator sameBins; if (sameBins.isValid(outputWs) != "") { throw std::invalid_argument("Axes must have common bins for this " "algorithm to work - try Rebin first"); } isRefAxis = true; } double axisValue(0); // Create muparser try { mu::Parser p; // set parameter lookups for the axis value, allow both cases p.DefineVar("y", &axisValue); p.DefineVar("x", &axisValue); p.DefineVar("Y", &axisValue); p.DefineVar("X", &axisValue); p.SetExpr(formula); try { if (isRefAxis) { int64_t numberOfSpectra_i = static_cast<int64_t>( outputWs->getNumberHistograms()); // cast to make openmp happy // Calculate the new (common) X values MantidVec::iterator iter; for (iter = outputWs->dataX(0).begin(); iter != outputWs->dataX(0).end(); ++iter) { axisValue = *iter; double result = p.Eval(); *iter = result; } auto xVals = outputWs->refX(0); Progress prog(this, 0.6, 1.0, numberOfSpectra_i); PARALLEL_FOR1(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) { // cppcheck cannot see that this is used by reference by muparser // cppcheck-suppress unreadVariable axisValue = axisPtr->getValue(i); double result = p.Eval(); axisPtr->setValue(i, result); } } } 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()); } } 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()); } IncreasingAxisValidator incAxis; std::string incAxisResult = incAxis.isValid(outputWs); if (incAxisResult != "") g_log.warning( incAxisResult + ".\n" "Some of the Mantid algorithms might not use the workspace correctly."); if ((axisUnits != "") || (axisTitle != "")) { if (axisTitle == "") { axisTitle = axisPtr->unit()->caption(); } if (axisUnits == "") { axisUnits = axisPtr->unit()->label(); } axisPtr->unit() = boost::make_shared<Units::Label>(axisTitle, axisUnits); } } } // namespace Algorithms } // namespace Mantid