//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include <cmath>
#include <sstream>
#include "MantidAlgorithms/BinaryOperation.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/WorkspaceIterator.h"
using namespace Mantid::API;
using namespace Mantid::Kernel;
namespace Mantid
{
namespace Algorithms
{
// Get a reference to the logger
Logger& BinaryOperation::g_log = Logger::get("BinaryOperation");
/** Initialisation method.
* Defines input and output workspaces
*
*/
void BinaryOperation::init()
{
declareProperty(new WorkspaceProperty<MatrixWorkspace>("InputWorkspace_1","",Direction::Input));
declareProperty(new WorkspaceProperty<MatrixWorkspace>("InputWorkspace_2","",Direction::Input));
declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace","",Direction::Output));
}
/** Executes the algorithm
*
* @throw runtime_error Thrown if algorithm cannot execute
*/
void BinaryOperation::exec()
{
// get input workspace, dynamic cast not needed
MatrixWorkspace_sptr in_work1 = getProperty("InputWorkspace_1");
MatrixWorkspace_sptr in_work2 = getProperty("InputWorkspace_2");
// Check that the input workspace are compatible
if (!checkCompatibility(in_work1,in_work2))
{
std::ostringstream ostr;
ostr << "The two workspaces are not compatible for algorithm " << this->name();
g_log.error() << ostr << std::endl;
throw std::invalid_argument( ostr.str() );
}
MatrixWorkspace::const_iterator ti_in1 = createConstIterator(in_work1,in_work2);
MatrixWorkspace::const_iterator ti_in2 = createConstIterator(in_work2,in_work1);
MatrixWorkspace_sptr out_work = getProperty("OutputWorkspace");
// 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 ( out_work == in_work1 )
{
if ( in_work2->size() > in_work1->size() ) out_work = createOutputWorkspace(in_work1,in_work2);
}
else if ( out_work == in_work2 )
{
if ( in_work1->size() > in_work2->size() ) out_work = createOutputWorkspace(in_work1,in_work2);
}
else
{
out_work = createOutputWorkspace(in_work1,in_work2);
}
MatrixWorkspace::iterator ti_out(*out_work);
//perform the operation through an abstract call
performBinaryOperation(ti_in1,ti_in2,ti_out);
// Assign it to the output workspace property
setProperty("OutputWorkspace",out_work);
return;
}
const bool BinaryOperation::checkCompatibility(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs) const
{
Unit_sptr lhs_unit = Unit_sptr();
Unit_sptr rhs_unit = Unit_sptr();
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();
}
// Check the workspaces have the same units and distribution flag
if ( lhs_unit != rhs_unit || lhs->YUnit() != rhs->YUnit() || lhs->isDistribution() != rhs->isDistribution() )
{
return false;
}
// Check the size compatibility
if (!checkSizeCompatibility(lhs,rhs))
{
std::ostringstream ostr;
ostr<<"The sizes of the two workspaces are not compatible for algorithm "<<this->name();
g_log.error() << ostr << std::endl;
throw std::invalid_argument( ostr.str() );
}
return true;
}
/** 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 true The two workspaces are size compatible
* @retval false The two workspaces are NOT size compatible
*/
const bool BinaryOperation::checkSizeCompatibility(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs) const
{
//in order to be size compatible then the larger workspace
//must divide by the size of the smaller workspace leaving no remainder
if (rhs->size() ==0) return false;
return ((lhs->size() % rhs->size()) == 0);
}
/** Performs a simple check to see if the X arrays of two workspaces are compatible for a binary operation
* The X arrays of two workspaces must be identical to allow a binary operation to be performed
* @param lhs the first workspace to compare
* @param rhs the second workspace to compare
* @retval true The two workspaces are size compatible
* @retval false The two workspaces are NOT size compatible
*/
const bool BinaryOperation::checkXarrayCompatibility(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs) const
{
// Not using the WorkspaceHelpers::matching bins method because that requires the workspaces to be
// the same size, which isn't a requirement of BinaryOperation
// single values, or workspaces with just a single bin/value in each spectrum, are compatible with anything
if ((rhs->blocksize() ==1) || (lhs->blocksize() ==1)) return true;
const std::vector<double>& w1x = lhs->readX(0);
const std::vector<double>& w2x = rhs->readX(0);
double sum;
sum=0.0;
for (unsigned int i=0; i < w1x.size(); i++) sum += fabs(w1x[i]-w2x[i]);
if( sum < 0.0000001)
return true;
else
return false;
}
/** Gets the number of time an iterator over the first workspace would have to loop to perform a full iteration of the second workspace
* @param lhs the first workspace to compare
* @param rhs the second workspace to compare
* @returns Integer division of rhs.size()/lhs.size() with a minimum of 1
*/
const int BinaryOperation::getRelativeLoopCount(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs) const
{
int lhsSize = lhs->size();
if (lhsSize == 0) return 1;
int retVal = rhs->size()/lhsSize;
return (retVal == 0)?1:retVal;
}
/** Creates a suitable output workspace for a binary operatiion based on the two input workspaces
* @param lhs the first workspace to compare
* @param rhs the second workspace to compare
* @returns a pointer to a new zero filled workspace the same type and size as the larger of the two input workspaces.
*/
API::MatrixWorkspace_sptr BinaryOperation::createOutputWorkspace(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs) const
{
//get the largest workspace
const API::MatrixWorkspace_const_sptr wsLarger = (lhs->size() > rhs->size()) ? lhs : rhs;
//create a new workspace
API::MatrixWorkspace_sptr retVal = API::WorkspaceFactory::Instance().create(wsLarger);
return retVal;
}
/** Creates a const iterator taking into account loop
* @param wsMain The workspace theat the iterator will be created for
* @param wsComparison The workspace to be used for axes comparisons
* @returns a const iterator to wsMain, with loop count and orientation set appropriately
*/
MatrixWorkspace::const_iterator BinaryOperation::createConstIterator(const API::MatrixWorkspace_const_sptr wsMain, const API::MatrixWorkspace_const_sptr wsComparison) const
{
//get loop count
unsigned int loopDirection = LoopOrientation::Vertical;
int loopCount = getRelativeLoopCount(wsMain,wsComparison);
if (loopCount > 1)
{
loopDirection = getLoopDirection(wsMain,wsComparison);
}
else
{
if (!checkXarrayCompatibility(wsMain,wsComparison))
{
g_log.error("The x arrays of the workspaces are not identical");
throw std::invalid_argument("The x arrays of the workspaces are not identical");
}
}
MatrixWorkspace::const_iterator it(*wsMain,loopCount,loopDirection);
return it;
}
/** Determines the required loop direction for a looping iterator
* @param wsMain The workspace theat the iterator will be created for
* @param wsComparison The workspace to be used for axes comparisons
* @returns An value describing the orientation of the 1D workspace to be looped
* @retval 0 Horizontal - The number and contents of the X axis bins match
* @retval 1 Vertical - The number of detector elements match
*/
unsigned int BinaryOperation::getLoopDirection(const API::MatrixWorkspace_const_sptr wsMain, const API::MatrixWorkspace_const_sptr wsComparison) const
{
unsigned int retVal = LoopOrientation::Horizontal;
//check if the vertical sizes match
int wsMainArraySize = wsMain->size(); //this must be a 1D array for this to work
int wsComparisonArraySize = wsComparison->size()/wsComparison->blocksize();
if (wsMainArraySize == wsComparisonArraySize)
{
retVal = LoopOrientation::Vertical;
}
//check if Horizontial looping matches in length
if (wsMain->blocksize() == wsComparison->blocksize())
{
//it does, now check if the X arrays are compatible
if (!checkXarrayCompatibility(wsMain,wsComparison))
{
if(retVal == LoopOrientation::Horizontal)
{
//this is a problem, the lengths only match Horizontally but the data does not match
g_log.error("The x arrays of the workspaces are not identical");
throw std::invalid_argument("The x arrays of the workspaces are not identical");
}
}
else
{
//all is good in the world
retVal = LoopOrientation::Horizontal;
}
}
return retVal;
}
}
}