diff --git a/Code/Mantid/Algorithms/inc/MantidAlgorithms/BinaryOperation.h b/Code/Mantid/Algorithms/inc/MantidAlgorithms/BinaryOperation.h index 5c0349d3dc3228193ed57d7a610cbda26da9d457..7ea658bbac92203625bb5717d834b73fab78b248 100644 --- a/Code/Mantid/Algorithms/inc/MantidAlgorithms/BinaryOperation.h +++ b/Code/Mantid/Algorithms/inc/MantidAlgorithms/BinaryOperation.h @@ -83,6 +83,11 @@ namespace Mantid /// Checks the overall size compatibility of two workspaces virtual bool checkSizeCompatibility(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs) const; + + /// Checks if the spectra at the given index of either input workspace is masked. If so then the output spectra has zeroed data + /// and is also masked. The function returns true if further processing is not required on the spectra. + virtual bool propagateSpectraMask(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs, + const int index, API::MatrixWorkspace_sptr out); /** Carries out the binary operation on a single spectrum. * @param lhsX The X values, made available if required. @@ -129,10 +134,14 @@ namespace Mantid void doSingleSpectrum(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs,API::MatrixWorkspace_sptr out); void doSingleColumn(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs,API::MatrixWorkspace_sptr out); void do2D(const API::MatrixWorkspace_const_sptr lhs,const API::MatrixWorkspace_const_sptr rhs,API::MatrixWorkspace_sptr out); - - void propagateMasking(const API::MatrixWorkspace_const_sptr rhs, API::MatrixWorkspace_sptr out); - - API::Progress* m_progress; ///< Progress reporting + void propagateBinMasks(const API::MatrixWorkspace_const_sptr rhs, API::MatrixWorkspace_sptr out); + /// Apply masking requested by propagateSpectraMasks. + void applyMaskingToOutput(API::MatrixWorkspace_sptr out); + + /// A store for accumulated spectra that should be masked in the output workspace + std::vector<int> m_indicesToMask; + /// Progress reporting + API::Progress* m_progress; }; } // namespace Algorithms diff --git a/Code/Mantid/Algorithms/src/BinaryOperation.cpp b/Code/Mantid/Algorithms/src/BinaryOperation.cpp index 9a0abfb6a1ad061e9205bac6bf1d83391373da97..3243884868885adeae0b959b974ad5eb6e773518 100644 --- a/Code/Mantid/Algorithms/src/BinaryOperation.cpp +++ b/Code/Mantid/Algorithms/src/BinaryOperation.cpp @@ -4,7 +4,9 @@ #include "MantidAlgorithms/BinaryOperation.h" #include "MantidAPI/WorkspaceProperty.h" #include "MantidAPI/FrameworkManager.h" +#include "MantidGeometry/IDetector.h" +using namespace Mantid::Geometry; using namespace Mantid::API; using namespace Mantid::Kernel; using namespace Mantid::DataObjects; @@ -13,7 +15,7 @@ namespace Mantid { namespace Algorithms { - BinaryOperation::BinaryOperation() : API::PairedGroupAlgorithm(), m_progress(NULL) {} + BinaryOperation::BinaryOperation() : API::PairedGroupAlgorithm(), m_indicesToMask(), m_progress(NULL) {} BinaryOperation::~BinaryOperation() { @@ -100,14 +102,16 @@ namespace Mantid } else if ( rhs->blocksize() == 1 ) // Single column on rhs { + m_indicesToMask.reserve(out_work->getNumberHistograms()); doSingleColumn(lhs,rhs,out_work); } else // The two are both 2D and should be the same size { + m_indicesToMask.reserve(out_work->getNumberHistograms()); do2D(lhs,rhs,out_work); } - // Call the virtual unit manipulation method + applyMaskingToOutput(out_work); setOutputUnits(lhs,rhs,out_work); // Assign the result to the output workspace property @@ -200,6 +204,46 @@ namespace Mantid return ( lhs->blocksize() == rhs->blocksize() && ( rhsSpec==1 || lhs->getNumberHistograms() == rhsSpec ) ); } + //-------------------------------------------------------------------------------------------- + /** + * Checks if the spectra at the given index of either input workspace is masked. If so then the output spectra has zeroed data + * and is also masked. + * @param lhs A pointer to the left-hand operand + * @param rhs A pointer to the right-hand operand + * @param index The workspace index to check + * @param out A pointer to the output workspace + * @returns True if further processing is not required on the spectra, false if the binary operation should be performed. + */ + bool BinaryOperation::propagateSpectraMask(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs, + const int index, API::MatrixWorkspace_sptr out) + { + bool continueOp(true); + IDetector_sptr det_lhs, det_rhs; + try + { + det_lhs = lhs->getDetector(index); + det_rhs = rhs->getDetector(index); + } + catch(std::runtime_error &) + { + } + if( (det_lhs && det_lhs->isMasked()) || ( det_rhs && det_rhs->isMasked()) ) + { + continueOp = false; + //Zero the output data and ensure that the output spectra is masked. The masking is done outside of this + //loop modiying the parameter map in a multithreaded loop requires too much locking + m_indicesToMask.push_back(index); + MantidVec & yValues = out->dataY(index); + MantidVec & eValues = out->dataE(index); + MantidVec::const_iterator yend = yValues.end(); + for( MantidVec::iterator yit(yValues.begin()), eit(eValues.begin()); yit != yend; ++yit, ++eit) + { + (*yit) = 0.0; + (*eit) = 0.0; + } + } + return continueOp; + } //-------------------------------------------------------------------------------------------- /** Called when the rhs operand is a single value. @@ -241,7 +285,7 @@ namespace Mantid 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); + propagateBinMasks(rhs,out); // Pull out the rhs spectrum const MantidVec& rhsY = rhs->readY(0); @@ -282,7 +326,10 @@ namespace Mantid const double rhsE = rhs->readE(i)[0]; out->setX(i,lhs->refX(i)); - performBinaryOperation(lhs->readX(i),lhs->readY(i),lhs->readE(i),rhsY,rhsE,out->dataY(i),out->dataE(i)); + if( propagateSpectraMask(lhs, rhs, i, out) ) + { + 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 } @@ -298,7 +345,7 @@ namespace Mantid 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); + propagateBinMasks(rhs,out); // Loop over the spectra calling the virtual function for each one const int numHists = lhs->getNumberHistograms(); PARALLEL_FOR3(lhs,rhs,out) @@ -306,7 +353,10 @@ namespace Mantid { PARALLEL_START_INTERUPT_REGION out->setX(i,lhs->refX(i)); - performBinaryOperation(lhs->readX(i),lhs->readY(i),lhs->readE(i),rhs->readY(i),rhs->readE(i),out->dataY(i),out->dataE(i)); + if( propagateSpectraMask(lhs, rhs, i, out) ) + { + 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 } @@ -318,7 +368,7 @@ namespace Mantid * @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) + void BinaryOperation::propagateBinMasks(const API::MatrixWorkspace_const_sptr rhs, API::MatrixWorkspace_sptr out) { const int outHists = out->getNumberHistograms(); const int rhsHists = rhs->getNumberHistograms(); @@ -331,10 +381,42 @@ namespace Mantid 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); + } } } } + /** + * Apply the requested masking to the output workspace + * @param out The workspace to mask + */ + void BinaryOperation::applyMaskingToOutput(API::MatrixWorkspace_sptr out) + { + int nindices = static_cast<int>(m_indicesToMask.size()); + ParameterMap &pmap = out->instrumentParameters(); + PARALLEL_FOR1(out) + for(int i = 0; i < nindices; ++i) + { + PARALLEL_START_INTERUPT_REGION + + try + { + IDetector_sptr det_out = out->getDetector(m_indicesToMask[i]); + PARALLEL_CRITICAL(BinaryOperation_masking) + { + pmap.addBool(det_out.get(), "masked", true); + } + } + catch(std::runtime_error &) + { + } + + PARALLEL_END_INTERUPT_REGION + } + PARALLEL_CHECK_INTERUPT_REGION + } + } // namespace Algorithms } // namespace Mantid diff --git a/Code/Mantid/Algorithms/test/BinaryOperationTest.h b/Code/Mantid/Algorithms/test/BinaryOperationTest.h index 6f9557334fecca1810cef902002bbd8dbb47e37d..a4e5207c2a8b02c65bc5686b1f3848af9dad15f5 100644 --- a/Code/Mantid/Algorithms/test/BinaryOperationTest.h +++ b/Code/Mantid/Algorithms/test/BinaryOperationTest.h @@ -22,6 +22,14 @@ public: BinaryOpHelper() : BinaryOperation() {}; /// Destructor virtual ~BinaryOpHelper() {}; + + /// function to return a name of the algorithm, must be overridden in all algorithms + virtual const std::string name() const { return "BinaryOpHelper"; } + /// function to return a version of the algorithm, must be overridden in all algorithms + virtual int version() const { return 1; } + /// function to return a category of the algorithm. A default implementation is provided + virtual const std::string category() const {return "Helper";} + const bool checkSizeCompatibility(const MatrixWorkspace_sptr ws1,const MatrixWorkspace_sptr ws2) const { return BinaryOperation::checkSizeCompatibility(ws1,ws2); @@ -102,6 +110,57 @@ public: TS_ASSERT(helper.checkSizeCompatibility(work_in1,work_inEvent2)); } + void testMaskedSpectraPropagation() + { + const int sizex = 10,sizey=20; + std::set<int> masking; + masking.insert(0); + masking.insert(2); + masking.insert(4); + + MatrixWorkspace_sptr work_in1 = WorkspaceCreationHelper::Create2DWorkspace123(sizex,sizey, 0, masking); + MatrixWorkspace_sptr work_in2 = WorkspaceCreationHelper::Create2DWorkspace154(sizex,sizey); + + BinaryOpHelper helper; + helper.initialize(); + helper.setProperty("LHSWorkspace", work_in1); + helper.setProperty("RHSWorkspace", work_in2); + const std::string outputSpace("test"); + helper.setPropertyValue("OutputWorkspace", outputSpace); + helper.setRethrows(true); + helper.execute(); + + + TS_ASSERT(helper.isExecuted()); + + MatrixWorkspace_sptr output = boost::dynamic_pointer_cast<MatrixWorkspace>(AnalysisDataService::Instance().retrieve(outputSpace)); + TS_ASSERT(output); + + for( int i = 0; i < sizey; ++i ) + { + IDetector_sptr det; + try + { + det = output->getDetector(i); + } + catch(Kernel::Exception::NotFoundError&) + { + } + + TS_ASSERT(det); + if( !det ) TS_FAIL("No detector found"); + if( masking.count(i) == 0 ) + { + TS_ASSERT_EQUALS(det->isMasked(), false); + } + else + { + TS_ASSERT_EQUALS(det->isMasked(), true); + } + } + + } + }; #endif /*BINARYOPERATIONTEST_H_*/ diff --git a/Code/Mantid/Algorithms/test/DivideTest.h b/Code/Mantid/Algorithms/test/DivideTest.h index d5663d2db65400b4381229e109cfc0bde9c1e5c1..76973d92d4e563a8a44a4986db3888e2f3488e9f 100644 --- a/Code/Mantid/Algorithms/test/DivideTest.h +++ b/Code/Mantid/Algorithms/test/DivideTest.h @@ -12,11 +12,13 @@ #include "MantidDataObjects/Workspace1D.h" #include "MantidAPI/WorkspaceProperty.h" #include "MantidAPI/WorkspaceOpOverloads.h" - +#include <limits> +#include "MantidGeometry/IDetector.h" using namespace Mantid::API; using namespace Mantid::Kernel; using namespace Mantid::Algorithms; using namespace Mantid::DataObjects; +using namespace Mantid::Geometry; class DivideTest : public CxxTest::TestSuite { @@ -49,6 +51,17 @@ public: perfomTest(work_in1,work_in2); } + void testDivideWithMaskedSpectraProducesZeroes() + { + doDivideWithMaskedTest(false); + } + + void testDivideWithMaskedSpectraProducesZeroesWhenReplacingInputWorkspace() + { + doDivideWithMaskedTest(true); + } + + void testExec1D2D() { int sizex = 10,sizey=20; @@ -222,6 +235,75 @@ private: TS_ASSERT_DELTA(err3, work_out1->readE(i/work_in1->blocksize())[i%work_in1->blocksize()], 0.0001); } + void doDivideWithMaskedTest(bool replaceInput) + { + const int sizex = 10,sizey=20; + std::set<int> masking; + masking.insert(0); + masking.insert(2); + masking.insert(7); + + MatrixWorkspace_sptr work_in1 = WorkspaceCreationHelper::Create2DWorkspace123(sizex,sizey, 0, masking); + MatrixWorkspace_sptr work_in2 = WorkspaceCreationHelper::Create2DWorkspace154(sizex,sizey, 0, masking); + const std::string lhs("work_in1"), rhs("work_in2"); + AnalysisDataService::Instance().add(lhs, work_in1); + AnalysisDataService::Instance().add(rhs, work_in2); + + //Zero some data to test mask propagation + for( int j = 0; j < sizex; ++j ) + { + work_in1->dataY(0)[j] = 0.0; + work_in1->dataY(2)[j] = 0.0; + work_in1->dataY(7)[j] = 0.0; + + work_in2->dataY(0)[j] = 0.0; + work_in2->dataY(2)[j] = 0.0; + work_in2->dataY(7)[j] = 0.0; + } + + Divide helper; + helper.initialize(); + helper.setPropertyValue("LHSWorkspace", lhs); + helper.setPropertyValue("RHSWorkspace", rhs); + std::string outputSpace; + if( replaceInput ) + { + outputSpace = lhs; + } + else + { + outputSpace = "lhsOverRhs"; + } + helper.setPropertyValue("OutputWorkspace", lhs); + helper.execute(); + + TS_ASSERT(helper.isExecuted()); + + MatrixWorkspace_sptr output = boost::dynamic_pointer_cast<MatrixWorkspace>(AnalysisDataService::Instance().retrieve("work_in1")); + TS_ASSERT(output); + + for( int i = 0; i < sizey; ++i ) + { + IDetector_sptr det = output->getDetector(i); + TS_ASSERT(det); + if( !det ) TS_FAIL("No detector found"); + if( masking.count(i) == 0 ) + { + TS_ASSERT_EQUALS(det->isMasked(), false); + } + else + { + TS_ASSERT_EQUALS(det->isMasked(), true); + double yValue = output->readY(i)[0]; + TS_ASSERT_EQUALS(yValue, yValue ); + TS_ASSERT_DIFFERS(yValue, std::numeric_limits<double>::infinity() ); + } + } + AnalysisDataService::Instance().remove(lhs); + AnalysisDataService::Instance().remove(rhs); + if( !replaceInput ) AnalysisDataService::Instance().remove(outputSpace); + } + }; diff --git a/Code/Mantid/Algorithms/test/MedianDetectorTestTest.h b/Code/Mantid/Algorithms/test/MedianDetectorTestTest.h index c4a3ed325f927ea4c837585cfa6d177384eb655a..d8f8529bdf59082711decacffcd1973e370ce21b 100644 --- a/Code/Mantid/Algorithms/test/MedianDetectorTestTest.h +++ b/Code/Mantid/Algorithms/test/MedianDetectorTestTest.h @@ -13,6 +13,8 @@ #include "MantidDataHandling/LoadInstrument.h" #include <boost/shared_ptr.hpp> #include <boost/lexical_cast.hpp> +#include <Poco/File.h> +#include <Poco/Path.h> #include <cmath> #include <iostream> #include <sstream> @@ -70,7 +72,7 @@ public: const int numberOfSpectra = outputMat->getNumberHistograms(); TS_ASSERT_EQUALS(numberOfSpectra, (int)Nhist); // the numbers below are threshold values that were found by trail and error running these tests - const int firstGoodSpec = 37; + const int firstGoodSpec = 36; const int lastGoodSpec = 95; for (int lHist = 1; lHist < firstGoodSpec; lHist++) { @@ -195,12 +197,13 @@ public: TS_ASSERT_EQUALS ( fileLine, correct ) } testFile.close(); - remove(m_OFileName.c_str()); + Poco::File(m_OFileName.c_str()).remove(); } MedianDetectorTestTest() : m_IWSName("MedianDetectorTestInput"), - m_OFileName("MedianDetectorTestTestFile.txt") + m_OFileName("") { + m_OFileName = Poco::Path().absolute().resolve("MedianDetectorTestTestFile.txt").toString(); using namespace Mantid; // Set up a small workspace for testing Workspace_sptr space = WorkspaceFactory::Instance().create("Workspace2D",Nhist,11,10);