BinaryOperation.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/BinaryOperation.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/FrameworkManager.h"

using namespace Mantid::API;
using namespace Mantid::Kernel;

namespace Mantid
{
  namespace Algorithms
  {
    BinaryOperation::BinaryOperation() : API::Algorithm(), 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));
      declareProperty(new WorkspaceProperty<MatrixWorkspace>(inputPropName2(),"",Direction::Input));
      declareProperty(new WorkspaceProperty<MatrixWorkspace>(outputPropName(),"",Direction::Output));
    }

    /** Executes the algorithm
     *
     *  @throw runtime_error Thrown if algorithm cannot execute
     */
    void BinaryOperation::exec()
    {
      // get input workspace, dynamic cast not needed
      MatrixWorkspace_const_sptr lhs = getProperty(inputPropName1());
      MatrixWorkspace_const_sptr rhs = getProperty(inputPropName2());

      // Check that the input workspace are compatible
      if (!checkCompatibility(lhs,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() );
      }

      // Make sure the left hand workspace is the larger (this is fine if we got past the compatibility checks)
      if ( rhs->size() > lhs->size() )
      {
        MatrixWorkspace_const_sptr smaller = lhs;
        lhs = rhs;
        rhs = smaller;
      }
      
      MatrixWorkspace_sptr out_work = getProperty(outputPropName());
      // 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 != lhs && out_work != rhs) || ( out_work == rhs && ( lhs->size() > rhs->size() ) ) )
      {
        out_work = WorkspaceFactory::Instance().create(lhs);
      }
      // Initialise the progress reporting object
      m_progress = new Progress(this,0.0,1.0,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
      if ( rhs->size() == 1 ) // Single value workspace
      {
        doSingleValue(lhs,rhs,out_work);
      }
      else if ( rhs->getNumberHistograms() == 1 ) // Single spectrum on rhs
      {
        doSingleSpectrum(lhs,rhs,out_work);
      }
      else if ( rhs->blocksize() == 1 ) // Single column on rhs
      {
        doSingleColumn(lhs,rhs,out_work);
      }
      else // The two are both 2D and should be the same size
      {
        do2D(lhs,rhs,out_work);
      }
      
      // Assign the result to the output workspace property
      setProperty(outputPropName(),out_work);

      return;
    }

    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->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
      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.str() << 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
     */
    bool BinaryOperation::checkSizeCompatibility(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs) const
    {
      const int lhsSize = lhs->size();
      const int rhsSize = rhs->size();
      // A SingleValueWorkspace matches anything
      if ( rhsSize == 1 ) return true;
      // The rhs must not be smaller than the lhs
      if ( lhsSize < rhsSize ) return false;
      // 
      // 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 true;
      // Past this point, we require the X arrays to match. Note this only checks the first spectrum
      if ( !WorkspaceHelpers::matchingBins(lhs,rhs,true) ) return false;
      
      const int rhsSpec = rhs->getNumberHistograms();
      return ( lhs->blocksize() == rhs->blocksize() && ( rhsSpec==1 || lhs->getNumberHistograms() == rhsSpec ) );
    }

    /** Called when the rhs operand is a single value. 
     *  Loops over the lhs workspace calling the abstract binary operation function. 
     *  @param lhs The workspace which is the left hand operand
     *  @param rhs The workspace which is the right hand operand
     *  @param out The result workspace
     */
    void BinaryOperation::doSingleValue(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs,API::MatrixWorkspace_sptr out)
    {
      // 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 = rhs->readY(0)[0];
      const double rhsE = rhs->readE(0)[0];
      
      // Now loop over the spectra of the left hand side calling the virtual function
      const int numHists = lhs->getNumberHistograms();
      PARALLEL_FOR3(lhs,rhs,out)
      for (int i = 0; i < numHists; ++i)
      {
        PARALLEL_START_INTERUPT_REGION
        out->dataX(i) = lhs->readX(i);
        performBinaryOperation(lhs->readX(i),lhs->readY(i),lhs->readE(i),rhsY,rhsE,out->dataY(i),out->dataE(i));
        m_progress->report();
        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION
   }
   
    /** Called when the rhs operand is a single spectrum. 
     *  Loops over the lhs workspace calling the abstract binary operation function. 
     *  @param lhs The workspace which is the left hand operand
     *  @param rhs The workspace which is the right hand operand
     *  @param out The result workspace
     */
    void BinaryOperation::doSingleSpectrum(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs,API::MatrixWorkspace_sptr out)
    {
      // Propagate any masking first or it could mess up the numbers
      propagateMasking(rhs,out);

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

      // Now loop over the spectra of the left hand side calling the virtual function
      const int numHists = lhs->getNumberHistograms();
      PARALLEL_FOR3(lhs,rhs,out)
      for (int i = 0; i < numHists; ++i)
      {
        PARALLEL_START_INTERUPT_REGION
        out->dataX(i) = lhs->readX(i);
        performBinaryOperation(lhs->readX(i),lhs->readY(i),lhs->readE(i),rhsY,rhsE,out->dataY(i),out->dataE(i));
        m_progress->report();
        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION
    }
    
    /** Called when the rhs operand is a 2D workspace of single values. 
     *  Loops over the workspaces calling the abstract binary operation function. 
     *  @param lhs The workspace which is the left hand operand
     *  @param rhs The workspace which is the right hand operand
     *  @param out The result workspace
     */
    void BinaryOperation::doSingleColumn(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs,API::MatrixWorkspace_sptr out)
    {
      // Don't propate masking from the rhs here - it would be decidedly odd if the single bin was masked

      // Now loop over the spectra of the left hand side pulling out the single value from each rhs 'spectrum'
      // and then calling the virtual function
      const int numHists = lhs->getNumberHistograms();
      PARALLEL_FOR3(lhs,rhs,out)
      for (int i = 0; i < numHists; ++i)
      {
        PARALLEL_START_INTERUPT_REGION
        const double rhsY = rhs->readY(i)[0];
        const double rhsE = rhs->readE(i)[0];        
        
        out->dataX(i) = lhs->readX(i);
        performBinaryOperation(lhs->readX(i),lhs->readY(i),lhs->readE(i),rhsY,rhsE,out->dataY(i),out->dataE(i));
        m_progress->report();
        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 lhs The workspace which is the left hand operand
     *  @param rhs The workspace which is the right hand operand
     *  @param out The result workspace
     */
    void BinaryOperation::do2D(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs,API::MatrixWorkspace_sptr out)
    {
      // Propagate any masking first or it could mess up the numbers
      propagateMasking(rhs,out);
      // Loop over the spectra calling the virtual function for each one
      const int numHists = lhs->getNumberHistograms();
      PARALLEL_FOR3(lhs,rhs,out)
      for (int i = 0; i < numHists; ++i)
      {
        PARALLEL_START_INTERUPT_REGION
        out->dataX(i) = lhs->readX(i);
        performBinaryOperation(lhs->readX(i),lhs->readY(i),lhs->readE(i),rhs->readY(i),rhs->readE(i),out->dataY(i),out->dataE(i));
        m_progress->report();
        PARALLEL_END_INTERUPT_REGION
      }      
      PARALLEL_CHECK_INTERUPT_REGION
    }

    /** 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::propagateMasking(const API::MatrixWorkspace_const_sptr rhs, API::MatrixWorkspace_sptr out)
    {
      const int outHists = out->getNumberHistograms();
      const int rhsHists = rhs->getNumberHistograms();
      for (int 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->maskBin(i,it->first,it->second);
        }
      }
    }

    /** This method is called if one of the selected workspaces for binary operation is a workspacegroup
     *  @param inputWSGrp pointer to the first workspace group
     *  @param prop a vector holding properties
     *  @retval false if  selected workspace groups sizes not match
     */
    bool BinaryOperation::processGroups(WorkspaceGroup_sptr inputWSGrp,const std::vector<Mantid::Kernel::Property*>&prop)
    {
      try
      {	
        //getting the input workspace group names
        const std::vector<std::string>& inputWSNames=inputWSGrp->getNames();
        int nSize=inputWSNames.size();
        //return if atleast one meber is not there in the group to process
        if(nSize<2)
        {
          g_log.error()<<"Input WorkspaceGroup has no members to process  "<<std::endl;
          return false;
        }
        std::vector<std::string> lhsWSGrp;
        std::vector<std::string> rhsWSGrp;
        bool bgroupExecStatus=true;
        //flag used for checking failure of  all memebrs of the group
        bool bgroupFailed=false;
        WorkspaceGroup_sptr outWSGrp= WorkspaceGroup_sptr(new WorkspaceGroup);

        getGroupNames(prop,lhsWSGrp,rhsWSGrp);
        //get member from  each group and call binary execute on each member
        IAlgorithm* alg = API::FrameworkManager::Instance().createAlgorithm(this->name() ,"",1);
        if(!alg)
        {
          g_log.error()<<"createAlgorithm failed for  "<<this->name()<<std::endl;
          return false;
        }
        std::vector<std::string>::const_iterator lhsItr=lhsWSGrp.begin(); 
        std::vector<std::string>::const_iterator rhsItr=rhsWSGrp.begin(); 
        int nPeriod=0;
        bool bStatus=0;

        if(lhsWSGrp.size()>1 && rhsWSGrp.size()>1)
        {
          if(isCompatibleSizes(lhsWSGrp,rhsWSGrp))
          {
            for (++lhsItr,++rhsItr;lhsItr!=lhsWSGrp.end();++lhsItr,++rhsItr)
            {
              ++nPeriod;
              setProperties(alg,prop,(*lhsItr),(*rhsItr),nPeriod,outWSGrp);
              bStatus=alg->execute();
              bgroupExecStatus=bgroupExecStatus&bStatus;
              bgroupFailed=bgroupFailed || bStatus;
            }
          }
        }
        else if(lhsWSGrp.size()==1)//if LHS is not a group workspace and RHS is group workspace
        {				
          for (++rhsItr;rhsItr!=rhsWSGrp.end();rhsItr++)
          {	++nPeriod;
          setProperties(alg,prop,(*lhsItr),(*rhsItr),nPeriod,outWSGrp);
          bStatus=alg->execute();
          bgroupExecStatus=bgroupExecStatus&bStatus;
          bgroupFailed=bgroupFailed || bStatus;
          }
        }
        else if (rhsWSGrp.size()==1)//if RHS is not a group workspace and LHS is a group workspace
        {				
          for (++lhsItr;lhsItr!=lhsWSGrp.end();lhsItr++)
          {	++nPeriod;
          setProperties(alg,prop,(*lhsItr),(*rhsItr),nPeriod,outWSGrp);
          bStatus=alg->execute();
          bgroupExecStatus=bgroupExecStatus && bStatus;
          bgroupFailed=bgroupFailed || bStatus;
          }
        }
        if(bgroupExecStatus)
          setExecuted(true);
        if(!bgroupFailed)
        {
          std::vector<std::string> names=outWSGrp->getNames();
          if(!names.empty())
            AnalysisDataService::Instance().remove(names[0]);
        }

        m_notificationCenter.postNotification(new FinishedNotification(this,this->isExecuted()));
        return bStatus;

      }
      catch(Exception::NotFoundError& )
      {
        return false;
      }
      catch(std::exception&)
      {
        return false;
      }
      return true;
    }

    /** This method sets properties for the algorithm
     *  @param alg pointer to the algorithm
     *  @param prop a vector holding properties
     *  @param lhsWSName name of the LHS workspace
     *  @param rhsWSName name of the RHS workspace
     *  @param nPeriod period number
     *  @param outWSGrp shared pointer to output workspace
     */
    void BinaryOperation::setProperties(IAlgorithm* alg,const std::vector<Kernel::Property*>&prop,
      const std::string& lhsWSName,const std::string& rhsWSName,int nPeriod,WorkspaceGroup_sptr outWSGrp)
    {
      std::string prevPropName("");
      std::string currentPropName("");

      std::vector<Mantid::Kernel::Property*>::const_iterator propItr=prop.begin();
      for (propItr=prop.begin();propItr!=prop.end();propItr++)
      {
        if(isWorkspaceProperty(*propItr))
        {
          if(isInputWorkspaceProperty(*propItr))
          {
            currentPropName=(*propItr)->name();
            if(prevPropName.empty())
            {	//LHS workspace 
              alg->setPropertyValue(currentPropName,lhsWSName);
            }
            else
            {	// RHS workspace 
              if(currentPropName.compare(prevPropName))
              { alg->setPropertyValue(currentPropName,rhsWSName);
              }
            }
            prevPropName=currentPropName;
          }//end of if loop for input workspace property
          if(isOutputWorkspaceProperty(*propItr))
          {
            std::string str=(*propItr)->value();
            std::stringstream speriodNum;
            speriodNum<<nPeriod;
            std::string outWSName=str+"_"+speriodNum.str();
            alg->setPropertyValue((*propItr)->name(),outWSName);
            if(nPeriod==1){
              if(outWSGrp)outWSGrp->add(str);
              AnalysisDataService::Instance().addOrReplace(str,outWSGrp );
            }
            if(outWSGrp)
              outWSGrp->add(outWSName);
          }//end of if loop for output workspace property
        }//end of if loop for workspace property
        else
        {
          alg->setPropertyValue((*propItr)->name(),(*propItr)->name());
        }
      }//end of for loop for property vector

    }

    /** This method checks both LHS and RHS workspaces are of same size
     *  @param lhsWSGrpNames a vector holding names of LHS Workspace
     *  @param rhsWSGrpNames a vector holding names of RHS Workspace
     *  @retval true if  selected workspace groups are of same size
     *  @retval false if  selected workspace groups sizes not match
     */
    bool BinaryOperation::isCompatibleSizes(const std::vector<std::string> &lhsWSGrpNames, const std::vector<std::string> &rhsWSGrpNames) const
    {
      //if both lhs and rhs workspaces are group workspaces ,check it's size
      if(lhsWSGrpNames.size()!=rhsWSGrpNames.size())
      {
        g_log.error("Selected workspace groups are not of same size to perform binary operation.");
        return false;
      }
      return true;
    }

    /** This method iterates through property vector and returns  LHS and RHS workspaces group names vectors
     *  @param prop  vector holding the properties
     *  @param lhsWSGrpNames a vector holding names of LHS Workspace
     *  @param rhsWSGrpNames a vector holding names of RHS Workspace
     */
    void BinaryOperation::getGroupNames(const std::vector<Kernel::Property*>&prop, 
      std::vector<std::string> &lhsWSGrpNames, std::vector<std::string> &rhsWSGrpNames) const
    {
      std::vector<Mantid::Kernel::Property*>::const_iterator itr;
      std::string prevPropName("");
      std::string currentPropName("");
      for (itr=prop.begin();itr!=prop.end();itr++)
      {	
        if(isWorkspaceProperty(*itr))
        {
          if(isInputWorkspaceProperty(*itr))
          {
            currentPropName=(*itr)->name();
            std::string wsName=(*itr)->value();
            Workspace_sptr wsPtr=AnalysisDataService::Instance().retrieve(wsName);
            WorkspaceGroup_sptr wsGrpSptr=boost::dynamic_pointer_cast<WorkspaceGroup>(wsPtr);
            if(prevPropName.empty())
            {	//LHSWorkspace 
              if(wsGrpSptr)
                lhsWSGrpNames=wsGrpSptr->getNames();
              else
                lhsWSGrpNames.push_back(wsName);
            }
            else
            {	
              if(currentPropName.compare(prevPropName))
              { //second workspace("RHSWorkSpace") 
                if(wsGrpSptr)
                  rhsWSGrpNames=wsGrpSptr->getNames();
                else
                  rhsWSGrpNames.push_back(wsName);
              }
            }
            prevPropName=currentPropName;
          }//end of if loop for input workspace property iteration
        }//end of if loop for workspace property iteration
      } //end of for loop for property iteration
    }

  } // namespace Algorithms
} // namespace Mantid