diff --git a/Code/Mantid/Framework/Algorithms/CMakeLists.txt b/Code/Mantid/Framework/Algorithms/CMakeLists.txt index d15239ba35c421b44c63a964d0e9646a20f4bd42..8e10e55664bac224bad533bdeac91427e7074b4a 100644 --- a/Code/Mantid/Framework/Algorithms/CMakeLists.txt +++ b/Code/Mantid/Framework/Algorithms/CMakeLists.txt @@ -191,6 +191,7 @@ set ( SRC_FILES src/UnwrapMonitor.cpp src/UnwrapSNS.cpp src/WeightedMean.cpp + src/WeightedMeanOfWorkspace.cpp src/WeightingStrategy.cpp src/XDataConverter.cpp ) @@ -389,6 +390,7 @@ set ( INC_FILES inc/MantidAlgorithms/UnwrapMonitor.h inc/MantidAlgorithms/UnwrapSNS.h inc/MantidAlgorithms/WeightedMean.h + inc/MantidAlgorithms/WeightedMeanOfWorkspace.h inc/MantidAlgorithms/WeightingStrategy.h inc/MantidAlgorithms/XDataConverter.h ) @@ -560,6 +562,7 @@ set ( TEST_FILES UnaryOperationTest.h UnwrapSNSTest.h UnwrapTest.h + WeightedMeanOfWorkspaceTest.h WeightedMeanTest.h WeightingStrategyTest.h WorkspaceCreationHelperTest.h diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/WeightedMeanOfWorkspace.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/WeightedMeanOfWorkspace.h new file mode 100644 index 0000000000000000000000000000000000000000..4e62cc3029c51db169aa09b81e7301e76cec0edf --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/WeightedMeanOfWorkspace.h @@ -0,0 +1,55 @@ +#ifndef MANTID_ALGORITHMS_WEIGHTEDMEANOFWORKSPACE_H_ +#define MANTID_ALGORITHMS_WEIGHTEDMEANOFWORKSPACE_H_ + +#include "MantidKernel/System.h" +#include "MantidAPI/Algorithm.h" + +namespace Mantid +{ + namespace Algorithms + { + + /** WeightedMeanOfWorkspace This algorithm calculates the weighted mean for + * a single workspace from all the detector (non-monitor, not masked) spectra + * in that workspace. The results is a single value for the entire workspace. + + Copyright © 2012 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. + + File change history is stored at: <https://svn.mantidproject.org/mantid/trunk/Code/Mantid> + Code Documentation is available at: <http://doxygen.mantidproject.org> + */ + class DLLExport WeightedMeanOfWorkspace : public API::Algorithm + { + public: + WeightedMeanOfWorkspace(); + virtual ~WeightedMeanOfWorkspace(); + + virtual const std::string name() const; + virtual int version() const; + virtual const std::string category() const; + + private: + virtual void initDocs(); + void init(); + void exec(); + }; + + } // namespace Algorithms +} // namespace Mantid + +#endif /* MANTID_ALGORITHMS_WEIGHTEDMEANOFWORKSPACE_H_ */ diff --git a/Code/Mantid/Framework/Algorithms/src/WeightedMeanOfWorkspace.cpp b/Code/Mantid/Framework/Algorithms/src/WeightedMeanOfWorkspace.cpp new file mode 100644 index 0000000000000000000000000000000000000000..196172785b174cc74915e62b0ba30cd30cad8677 --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/src/WeightedMeanOfWorkspace.cpp @@ -0,0 +1,131 @@ +/*WIKI* + +This algorithm calculates the weighted mean from all the spectra in a given +workspace. Monitors and masked spectra are ignored. Also, individual bins with +IEEE values will be excluded from the result. The weighted mean calculated by +the following: + +<math>\displaystyle y=\frac{\sum\frac{x_i}{\sigma^{2}_i}}{\sum\frac{1}{\sigma^{2}_i}} </math> + +and the variance is calculated by: + +<math>\displaystyle \sigma^{2}_y=\frac{1}{\sum\frac{1}{\sigma^{2}_i}} </math> + +*WIKI*/ + +#include "MantidAlgorithms/WeightedMeanOfWorkspace.h" +#include "MantidDataObjects/EventWorkspace.h" + +namespace Mantid +{ + using namespace API; + using namespace DataObjects; + using namespace Geometry; + using namespace Kernel; + + namespace Algorithms + { + // Register the algorithm into the AlgorithmFactory + DECLARE_ALGORITHM(WeightedMeanOfWorkspace) + + //---------------------------------------------------------------------------------------------- + /** Constructor + */ + WeightedMeanOfWorkspace::WeightedMeanOfWorkspace() + { + } + + //---------------------------------------------------------------------------------------------- + /** Destructor + */ + WeightedMeanOfWorkspace::~WeightedMeanOfWorkspace() + { + } + + //---------------------------------------------------------------------------------------------- + /// Algorithm's name for identification. @see Algorithm::name + const std::string WeightedMeanOfWorkspace::name() const { return "WeightedMeanOfWorkspace"; }; + + /// Algorithm's version for identification. @see Algorithm::version + int WeightedMeanOfWorkspace::version() const { return 1; }; + + /// Algorithm's category for identification. @see Algorithm::category + const std::string WeightedMeanOfWorkspace::category() const { return "Arithmetic"; } + + //---------------------------------------------------------------------------------------------- + /// Sets documentation strings for this algorithm + void WeightedMeanOfWorkspace::initDocs() + { + this->setWikiSummary("This algorithm calculates the weighted mean for an entire workspace."); + this->setOptionalMessage("This algorithm calculates the weighted mean for an entire workspace."); + } + + //---------------------------------------------------------------------------------------------- + /** Initialize the algorithm's properties. + */ + void WeightedMeanOfWorkspace::init() + { + this->declareProperty(new WorkspaceProperty<>("InputWorkspace", "", + Direction::Input), "An input workspace."); + this->declareProperty(new WorkspaceProperty<>("OutputWorkspace", "", + Direction::Output), "An output workspace."); + } + + //---------------------------------------------------------------------------------------------- + /** Execute the algorithm. + */ + void WeightedMeanOfWorkspace::exec() + { + MatrixWorkspace_sptr inputWS = this->getProperty("InputWorkspace"); + g_log.warning() << "Help1" << std::endl; + // Check if it is an event workspace + EventWorkspace_const_sptr eventW = boost::dynamic_pointer_cast<const EventWorkspace>(inputWS); + if (eventW != NULL) + { + throw std::runtime_error("WeightedMeanOfWorkspace cannot handle EventWorkspaces!"); + } + g_log.warning() << "Help2" << std::endl; + // Create the output workspace + MatrixWorkspace_sptr singleValued = WorkspaceFactory::Instance().create("WorkspaceSingleValue", 1, 1, 1); + g_log.warning() << "Help3" << std::endl; + // Calculate weighted mean + std::size_t numHists = inputWS->getNumberHistograms(); + double averageValue = 0.0; + double weightSum = 0.0; + for (std::size_t i = 0; i < numHists; ++i) + { + g_log.warning() << "Iteration " << i << std::endl; + IDetector_const_sptr det = inputWS->getDetector(i); + if( det->isMonitor() || det->isMasked() ) + { + continue; + } + g_log.warning() << "OK Pixel" << std::endl; + MantidVec y = inputWS->dataY(i); + MantidVec e = inputWS->dataE(i); + double weight = 0.0; + g_log.warning() << "A: " << y.size() << std::endl; + for (std::size_t j = 0; j < y.size(); ++j) + { + if (!boost::math::isnan(y[j]) && !boost::math::isinf(y[j]) && + !boost::math::isnan(e[j]) && !boost::math::isinf(e[j])) + { + weight = 1.0 / (e[j] * e[j]); + g_log.warning() << "B: " << weight << std::endl; + averageValue += (y[j] * weight); + g_log.warning() << "G: " << averageValue << std::endl; + weightSum += weight; + } + } + } + g_log.warning() << "Help4 " << averageValue << ", " << weightSum << std::endl; + singleValued->dataX(0)[0] = 0.0; + singleValued->dataY(0)[0] = averageValue / weightSum; + singleValued->dataE(0)[0] = std::sqrt(weightSum); + g_log.warning() << "Help5" << std::endl; + this->setProperty("OutputWorkspace", singleValued); + g_log.warning() << "Help6" << std::endl; + } + + } // namespace Algorithms +} // namespace Mantid diff --git a/Code/Mantid/Framework/Algorithms/test/WeightedMeanOfWorkspaceTest.h b/Code/Mantid/Framework/Algorithms/test/WeightedMeanOfWorkspaceTest.h new file mode 100644 index 0000000000000000000000000000000000000000..131732da841215840a70601aef8f7127e0954e41 --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/test/WeightedMeanOfWorkspaceTest.h @@ -0,0 +1,126 @@ +#ifndef WEIGHTEDMEANOFWORKSPACETEST_H_ +#define WEIGHTEDMEANOFWORKSPACETEST_H_ + +#include <cxxtest/TestSuite.h> + +#include "MantidAlgorithms/WeightedMeanOfWorkspace.h" +#include "MantidAPI/AnalysisDataService.h" +#include "MantidDataObjects/EventWorkspace.h" +#include "MantidTestHelpers/WorkspaceCreationHelper.h" + +using namespace Mantid::API; +using namespace Mantid::Algorithms; +using namespace Mantid::DataObjects; + +class WeightedMeanOfWorkspaceTest : public CxxTest::TestSuite +{ +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static WeightedMeanOfWorkspaceTest *createSuite() { return new WeightedMeanOfWorkspaceTest(); } + static void destroySuite( WeightedMeanOfWorkspaceTest *suite ) { delete suite; } + + void testInit() + { + WeightedMeanOfWorkspace alg; + TS_ASSERT_THROWS_NOTHING( alg.initialize() ) + TS_ASSERT( alg.isInitialized() ) + } + + void testExec() + { + // Get input workspace + MatrixWorkspace_sptr inputWS = createWorkspace(); + // Name of the output workspace. + std::string outWSName("WeightedMeanOfWorkspaceTest_OutputWS"); + + WeightedMeanOfWorkspace alg; + TS_ASSERT_THROWS_NOTHING( alg.initialize() ) + TS_ASSERT( alg.isInitialized() ) + TS_ASSERT_THROWS_NOTHING( alg.setProperty("InputWorkspace", inputWS) ); + TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) ); + TS_ASSERT_THROWS_NOTHING( alg.execute(); ); + TS_ASSERT( alg.isExecuted() ); + + // Retrieve the workspace from data service. + MatrixWorkspace_sptr ws; + TS_ASSERT_THROWS_NOTHING( ws = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(outWSName) ); + TS_ASSERT(ws); + if (!ws) return; + + TS_ASSERT_EQUALS( ws->getNumberHistograms(), 1 ); + TS_ASSERT_EQUALS( ws->readY(0)[0], 2.0 ); + TS_ASSERT_EQUALS( ws->readE(0)[0], 1.0 ); + + // Remove workspace from the data service. + AnalysisDataService::Instance().remove(outWSName); + } + + void testBadValues() + { + // Get input workspace + MatrixWorkspace_sptr inputWS = createWorkspace(false); + + // Put bad values into workspace + inputWS->dataY(1)[0] = 0.0/0.0; + inputWS->dataE(1)[1] = 0.0/0.0; + inputWS->dataY(1)[2] = 1.0/0.0; + + // Name of the output workspace. + std::string outWSName("WeightedMeanOfWorkspaceTest_OutputWS"); + + WeightedMeanOfWorkspace alg; + alg.initialize(); + alg.isInitialized(); + alg.setProperty("InputWorkspace", inputWS); + alg.setPropertyValue("OutputWorkspace", outWSName); + alg.execute(); + alg.isExecuted(); + + // Retrieve the workspace from data service. + MatrixWorkspace_sptr ws; + ws = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(outWSName); + if (!ws) return; + + TS_ASSERT_EQUALS( ws->getNumberHistograms(), 1 ); + TS_ASSERT_EQUALS( ws->readY(0)[0], 2.0 ); + TS_ASSERT_EQUALS( ws->readE(0)[0], 1.0 ); + + // Remove workspace from the data service. + AnalysisDataService::Instance().remove(outWSName); + } + + void testEventWs() + { + // Get input workspace + EventWorkspace_sptr inputWS = createEventWorkspace(); + // Name of the output workspace. + std::string outWSName("WeightedMeanOfWorkspaceTest_OutputWS"); + + WeightedMeanOfWorkspace alg; + alg.initialize(); + alg.isInitialized(); + alg.setProperty("InputWorkspace", inputWS); + alg.setPropertyValue("OutputWorkspace", outWSName); + alg.execute(); + TS_ASSERT( !alg.isExecuted() ); + } + +private: + MatrixWorkspace_sptr createWorkspace(bool doMasked = true) + { + std::set<int64_t> masked; + if (doMasked) + { + masked.insert(0); + } + return WorkspaceCreationHelper::Create2DWorkspace123(4, 3, true, masked); + } + + EventWorkspace_sptr createEventWorkspace() + { + return WorkspaceCreationHelper::CreateEventWorkspace(); + } +}; + +#endif /* WEIGHTEDMEANOFWORKSPACETEST_H_ */