Skip to content

Commit

Permalink
refs #5626. Test Parallel and Sequential implementations.
Browse files Browse the repository at this point in the history
  • Loading branch information
@OwenArnold
OwenArnold committed Jul 19, 2012
1 parent 9a0e9ef commit a213408
Show file tree
Hide file tree
Showing 3 changed files with 162 additions and 73 deletions.
11 changes: 8 additions & 3 deletions Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/NormaliseByDetector.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,14 @@

namespace Mantid
{

namespace API
{
/// Forward declaration for MatrixWorkspace.
class MatrixWorkspace;
}
namespace Algorithms
{

/** NormaliseByDetector : Normalises a workspace with respect to the detector efficiency function stored against components in the instrument parameters. See wiki for more details.
Detector efficiency functions are calculated using the wavelengths in the input workspace.
Expand Down Expand Up @@ -44,18 +45,22 @@ namespace Algorithms
class DLLExport NormaliseByDetector : public API::Algorithm
{
public:
NormaliseByDetector();
NormaliseByDetector(bool parallelExecution = true);
virtual ~NormaliseByDetector();

virtual const std::string name() const;
virtual int version() const;
virtual const std::string category() const;

private:
/// Flag to indicate that the histograms should be processed in parallel.
const bool m_parallelExecution;
/// Try to parse a function parameter and extract the correctly typed parameter.
const Mantid::Geometry::FitParameter tryParseFunctionParameter(Mantid::Geometry::Parameter_sptr parameter, Geometry::IDetector_const_sptr det);
/// Block to process histograms.
void processHistograms(boost::shared_ptr<Mantid::API::MatrixWorkspace> denominatorWS, boost::shared_ptr<const Mantid::API::MatrixWorkspace> inWS);
boost::shared_ptr<Mantid::API::MatrixWorkspace> processHistograms(boost::shared_ptr<Mantid::API::MatrixWorkspace> inWS);
/// Process indivdual histogram.
void processHistogram(size_t wsIndex, boost::shared_ptr<Mantid::API::MatrixWorkspace> denominatorWS, boost::shared_ptr<const Mantid::API::MatrixWorkspace> inWS);
virtual void initDocs();
void init();
void exec();
Expand Down
154 changes: 90 additions & 64 deletions Code/Mantid/Framework/Algorithms/src/NormaliseByDetector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,12 +29,10 @@ namespace Mantid
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(NormaliseByDetector)



//----------------------------------------------------------------------------------------------
/** Constructor
*/
NormaliseByDetector::NormaliseByDetector()
NormaliseByDetector::NormaliseByDetector(bool parallelExecution) : m_parallelExecution(parallelExecution)
{
}

Expand Down Expand Up @@ -91,75 +89,109 @@ namespace Mantid
return parameter->value<Geometry::FitParameter>();
}


/**
Process each histogram of the input workspace, extracting the detector/component and looking up the efficiency function.
Efficiency functions are then executed against the X data of the input workspace to generate new Y and E outputs for the denominatorWS.
@param wsIndex: The index of the histogram in the input workspace to process.
@param denominatorWS : Workspace that will become the denominator in the normalisation routine.
@param inputWorkspace: Workspace input. Contains instrument to use as well as X data to use.
*/
void NormaliseByDetector::processHistograms(MatrixWorkspace_sptr denominatorWS, MatrixWorkspace_const_sptr inWS)
void NormaliseByDetector::processHistogram(size_t wsIndex, MatrixWorkspace_sptr denominatorWS, MatrixWorkspace_const_sptr inWS)
{
const Geometry::ParameterMap& paramMap = inWS->instrumentParameters();
PARALLEL_FOR2(inWS, denominatorWS)
for(int wsIndex = 0; wsIndex < inWS->getNumberHistograms(); ++wsIndex)
{
PARALLEL_START_INTERUPT_REGION
Geometry::IDetector_const_sptr det = inWS->getDetector( wsIndex );
const std::string type = "fitting";
Geometry::Parameter_sptr foundParam = paramMap.getRecursiveByType(&(*det), type);
Geometry::IDetector_const_sptr det = inWS->getDetector( wsIndex );
const std::string type = "fitting";
Geometry::Parameter_sptr foundParam = paramMap.getRecursiveByType(&(*det), type);

const Geometry::FitParameter& foundFittingParam = tryParseFunctionParameter(foundParam, det);
const Geometry::FitParameter& foundFittingParam = tryParseFunctionParameter(foundParam, det);

std::string fitFunctionName = foundFittingParam.getFunction();
IFunction_sptr function = FunctionFactory::Instance().createFunction(fitFunctionName);
typedef std::vector<std::string> ParamNames;
ParamNames allParamNames = function->getParameterNames();
std::string fitFunctionName = foundFittingParam.getFunction();
IFunction_sptr function = FunctionFactory::Instance().createFunction(fitFunctionName);
typedef std::vector<std::string> ParamNames;
ParamNames allParamNames = function->getParameterNames();

// Lookup each parameter name.
for(ParamNames::iterator it = allParamNames.begin(); it != allParamNames.end(); ++it)
{
Geometry::Parameter_sptr param = paramMap.getRecursive(&(*det), (*it), type);

const Geometry::FitParameter& fitParam = tryParseFunctionParameter(param, det);

if ( fitParam.getFormula().compare("") == 0 )
{
throw std::invalid_argument("A Forumla has not been provided for a fit function");
}
else
{
std::string resultUnitStr = fitParam.getResultUnit();
if ( !resultUnitStr.empty() && resultUnitStr.compare("Wavelength") != 0)
{
throw std::invalid_argument("Units for function parameters must be in Wavelength");
}
}
mu::Parser p;
p.SetExpr(fitParam.getFormula());
double paramValue = p.Eval();
//Set the function coeffiecents.
function->setParameter(fitParam.getName(), paramValue);
}
// Lookup each parameter name.
for(ParamNames::iterator it = allParamNames.begin(); it != allParamNames.end(); ++it)
{
Geometry::Parameter_sptr param = paramMap.getRecursive(&(*det), (*it), type);

const Geometry::FitParameter& fitParam = tryParseFunctionParameter(param, det);

auto wavelengths = inWS->readX(wsIndex);
const size_t nInputBins = wavelengths.size() -1;
std::vector<double> centerPointWavelength(nInputBins);
std::vector<double> outIntensity(nInputBins);
for(size_t binIndex = 0; binIndex < nInputBins; ++binIndex)
if ( fitParam.getFormula().compare("") == 0 )
{
throw std::invalid_argument("A Forumla has not been provided for a fit function");
}
else
{
std::string resultUnitStr = fitParam.getResultUnit();
if ( !resultUnitStr.empty() && resultUnitStr.compare("Wavelength") != 0)
{
centerPointWavelength[binIndex] = 0.5*(wavelengths[binIndex] + wavelengths[binIndex+1]);
}
FunctionDomain1DVector domain(centerPointWavelength);
FunctionValues values(domain);
function->function(domain,values);
for(size_t i = 0; i < domain.size(); ++i)
throw std::invalid_argument("Units for function parameters must be in Wavelength");
}
}
mu::Parser p;
p.SetExpr(fitParam.getFormula());
double paramValue = p.Eval();
//Set the function coeffiecents.
function->setParameter(fitParam.getName(), paramValue);
}

auto wavelengths = inWS->readX(wsIndex);
const size_t nInputBins = wavelengths.size() -1;
std::vector<double> centerPointWavelength(nInputBins);
std::vector<double> outIntensity(nInputBins);
for(size_t binIndex = 0; binIndex < nInputBins; ++binIndex)
{
centerPointWavelength[binIndex] = 0.5*(wavelengths[binIndex] + wavelengths[binIndex+1]);
}
FunctionDomain1DVector domain(centerPointWavelength);
FunctionValues values(domain);
function->function(domain,values);
for(size_t i = 0; i < domain.size(); ++i)
{
outIntensity[i] = values[i];
}
denominatorWS->dataY(wsIndex) = outIntensity;
}

/**
Controlling function. Processes the histograms either in parallel or sequentially.
@param denominatorWS : Workspace that will become the denominator in the normalisation routine.
@param inputWorkspace: Workspace input. Contains instrument to use as well as X data to use.
*/
MatrixWorkspace_sptr NormaliseByDetector::processHistograms(MatrixWorkspace_sptr inWS)
{
// Clone the input workspace to create a template for the denominator workspace.
IAlgorithm_sptr cloneAlg = this->createSubAlgorithm("CloneWorkspace", 0.0, 1.0, true);
cloneAlg->setProperty("InputWorkspace", inWS);
cloneAlg->setPropertyValue("OutputWorkspace", "temp");
cloneAlg->executeAsSubAlg();
Workspace_sptr temp = cloneAlg->getProperty("OutputWorkspace");
MatrixWorkspace_sptr denominatorWS = boost::dynamic_pointer_cast<MatrixWorkspace>(temp);

// Choose between parallel execution and sequential execution then, process histograms accordingly.
const size_t nHistograms = inWS->getNumberHistograms();
if(m_parallelExecution == true)
{
PARALLEL_FOR2(inWS, denominatorWS)
for(int wsIndex = 0; wsIndex < nHistograms; ++wsIndex)
{
outIntensity[i] = values[i];
PARALLEL_START_INTERUPT_REGION
this->processHistogram(wsIndex, denominatorWS, inWS);
PARALLEL_END_INTERUPT_REGION
}
denominatorWS->dataY(wsIndex) = outIntensity;
PARALLEL_END_INTERUPT_REGION
PARALLEL_CHECK_INTERUPT_REGION
}
else
{
for(size_t wsIndex = 0; wsIndex < nHistograms; ++wsIndex)
{
this->processHistogram(wsIndex, denominatorWS, inWS);
}
PARALLEL_CHECK_INTERUPT_REGION
}

return denominatorWS;
};

//----------------------------------------------------------------------------------------------
Expand All @@ -169,14 +201,8 @@ namespace Mantid
{
MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");

IAlgorithm_sptr cloneAlg = this->createSubAlgorithm("CloneWorkspace", 0.0, 1.0, true);
cloneAlg->setProperty("InputWorkspace", inWS);
cloneAlg->setPropertyValue("OutputWorkspace", "temp");
cloneAlg->executeAsSubAlg();
Workspace_sptr temp = cloneAlg->getProperty("OutputWorkspace");
MatrixWorkspace_sptr denominatorWS = boost::dynamic_pointer_cast<MatrixWorkspace>(temp);

processHistograms(denominatorWS, inWS);
// Do the work of extracting functions and applying them to each bin on each histogram. The denominator workspace is mutable.
MatrixWorkspace_sptr denominatorWS = processHistograms(inWS);

// Perform the normalisation.
IAlgorithm_sptr divideAlg = this->createSubAlgorithm("Divide", 0.0, 1.0, true);
Expand Down
70 changes: 64 additions & 6 deletions Code/Mantid/Framework/Algorithms/test/NormaliseByDetectorTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -253,9 +253,9 @@ class NormaliseByDetectorTest : public CxxTest::TestSuite
/**
Helper method for running the algorithm and testing for invalid argument on execution.
*/
void do_test_throws_invalid_argument_on_execution(MatrixWorkspace_sptr inputWS)
void do_test_throws_invalid_argument_on_execution(MatrixWorkspace_sptr inputWS, bool parallel = true)
{
NormaliseByDetector alg;
NormaliseByDetector alg(parallel);
alg.setRethrows(true);
alg.initialize();
alg.setPropertyValue("OutputWorkspace", "out");
Expand All @@ -266,9 +266,9 @@ class NormaliseByDetectorTest : public CxxTest::TestSuite
/**
Helper method for running the algorithm and simply verifying that it runs without exception producing an output workspace..
*/
MatrixWorkspace_sptr do_test_doesnt_throw_on_execution(MatrixWorkspace_sptr inputWS)
MatrixWorkspace_sptr do_test_doesnt_throw_on_execution(MatrixWorkspace_sptr inputWS, bool parallel = true)
{
NormaliseByDetector alg;
NormaliseByDetector alg(parallel);
alg.setRethrows(true);
alg.initialize();
alg.setPropertyValue("OutputWorkspace", "out");
Expand Down Expand Up @@ -322,11 +322,20 @@ class NormaliseByDetectorTest : public CxxTest::TestSuite
do_test_doesnt_throw_on_execution(inputWS2);
}

void test_applies_instrument_function_to_child_detectors_throws_nothing()
void test_parallel_application_throws_nothing()
{
// Linear function 2*x + 1 applied to each x-value.
MatrixWorkspace_sptr inputWS = create_workspace_with_fitting_functions();
const bool parallel = true;
do_test_doesnt_throw_on_execution(inputWS, parallel);
}

void test_sequential_application_throws_nothing()
{
// Linear function 2*x + 1 applied to each x-value.
MatrixWorkspace_sptr inputWS = create_workspace_with_fitting_functions();
do_test_doesnt_throw_on_execution(inputWS);
const bool parallel = false;
do_test_doesnt_throw_on_execution(inputWS, parallel);
}

void test_workspace_with_instrument_only_fitting_functions()
Expand Down Expand Up @@ -361,6 +370,55 @@ class NormaliseByDetectorTest : public CxxTest::TestSuite
}
}

void test_compare_sequential_and_parallel_results()
{
const std::string outWSName = "normalised_ws";
// Linear function 2*x + 1 applied to each x-value. INSTRUMENT LEVEL FIT FUNCTION ONLY.
MatrixWorkspace_sptr inputWS = create_workspace_with_fitting_functions();
// Extract the output workspace so that we can verify the normalisation.
const bool parallel = true;
MatrixWorkspace_sptr outWS_parallel = do_test_doesnt_throw_on_execution(inputWS, parallel); //EXECUTES THE ALG IN PARALLEL.
MatrixWorkspace_sptr outWS_sequential = do_test_doesnt_throw_on_execution(inputWS, !parallel); //EXECUTES THE ALG SEQUENTIALLY.

// Output workspaces should have same number of histograms.
TS_ASSERT_EQUALS(2, outWS_parallel->getNumberHistograms());
TS_ASSERT_EQUALS(outWS_parallel->getNumberHistograms(), outWS_sequential->getNumberHistograms());

// Test the application of the linear function
for(size_t wsIndex = 0; wsIndex < inputWS->getNumberHistograms(); ++wsIndex)
{
const MantidVec& yValuesParallel = outWS_parallel->readY(wsIndex);
const MantidVec& xValuesParallel = outWS_parallel->readX(wsIndex);
const MantidVec& eValuesParallel = outWS_parallel->readE(wsIndex);

const MantidVec& yValuesSequential = outWS_sequential->readY(wsIndex);
const MantidVec& xValuesSequential = outWS_sequential->readX(wsIndex);
const MantidVec& eValuesSequential = outWS_sequential->readE(wsIndex);

// Compare against known sizes.
TS_ASSERT_EQUALS(3, yValuesParallel.size());
TS_ASSERT_EQUALS(3, eValuesParallel.size());
TS_ASSERT_EQUALS(4, xValuesParallel.size());
// Compare results from different execution types.
TS_ASSERT_EQUALS(yValuesSequential.size(), yValuesParallel.size());
TS_ASSERT_EQUALS(xValuesSequential.size(), xValuesParallel.size());
TS_ASSERT_EQUALS(eValuesSequential.size(), eValuesParallel.size());

const MantidVec& yInputValues = inputWS->readY(wsIndex);
const MantidVec& xInputValues = inputWS->readX(wsIndex);

for(size_t binIndex = 0; binIndex < (xInputValues.size() - 1); ++binIndex)
{
const double wavelength = (xInputValues[binIndex] + xInputValues[binIndex+1])/2;
const double expectedValue = yInputValues[binIndex] / ( (2*wavelength) + 1 ); // According to the equation written into the instrument parameter file for the instrument component link.
// Compare against the known/calculated value.
TS_ASSERT_EQUALS(expectedValue, yValuesParallel[binIndex]);
// Compare results from different execution types.
TS_ASSERT_EQUALS(yValuesSequential[binIndex], yValuesParallel[binIndex]);
}
}
}

void test_workspace_with_detector_level_only_fit_functions()
{
const std::string outWSName = "normalised_ws";
Expand Down

0 comments on commit a213408

Please sign in to comment.