Skip to content
Snippets Groups Projects
Commit 025f5b57 authored by Michael Reuter's avatar Michael Reuter
Browse files

Refs #5890. First cut at algorithm.

parent 30aaea27
No related merge requests found
......@@ -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
......
#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_ */
/*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
#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_ */
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment