Skip to content

Commit

Permalink
refs #11393. Normalization workspace option.
Browse files Browse the repository at this point in the history
  • Loading branch information
@OwenArnold
OwenArnold committed Mar 31, 2015
1 parent c66bd0c commit ca69654
Show file tree
Hide file tree
Showing 4 changed files with 355 additions and 153 deletions.
3 changes: 2 additions & 1 deletion Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/SmoothMD.h
View file
Expand Up @@ -4,6 +4,7 @@
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
#include <boost/shared_ptr.hpp>
#include <boost/optional.hpp>

namespace Mantid
{
Expand Down Expand Up @@ -49,7 +50,7 @@ namespace MDAlgorithms
std::map<std::string, std::string> validateInputs();

boost::shared_ptr<Mantid::API::IMDHistoWorkspace> hatSmooth(boost::shared_ptr<const Mantid::API::IMDHistoWorkspace> toSmooth,
const std::vector<int> &widthVector);
const std::vector<int> &widthVector, boost::optional<boost::shared_ptr<const Mantid::API::IMDHistoWorkspace> > weighting);

private:

Expand Down
100 changes: 92 additions & 8 deletions Code/Mantid/Framework/MDAlgorithms/src/SmoothMD.cpp
View file
Expand Up @@ -17,6 +17,7 @@
#include <string>
#include <sstream>
#include <utility>
#include <limits>
#include <boost/function.hpp>
#include <boost/bind.hpp>
#include <boost/scoped_ptr.hpp>
Expand All @@ -26,9 +27,17 @@ using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::MDEvents;

// Typedef for with vector
typedef std::vector<int> WidthVector;

// Typedef for an optional md histo workspace
typedef boost::optional<IMDHistoWorkspace_const_sptr> OptionalIMDHistoWorkspace_const_sptr;

// Typedef for a smoothing function
typedef boost::function<IMDHistoWorkspace_sptr(
IMDHistoWorkspace_const_sptr, const WidthVector &)> SmoothFunction;
IMDHistoWorkspace_const_sptr, const WidthVector &, OptionalIMDHistoWorkspace_const_sptr)> SmoothFunction;

// Typedef for a smoothing function map keyed by name.
typedef std::map<std::string, SmoothFunction> SmoothFunctionMap;

namespace {
Expand All @@ -50,7 +59,7 @@ std::vector<std::string> functions() {
*/
SmoothFunctionMap makeFunctionMap(Mantid::MDAlgorithms::SmoothMD * instance) {
SmoothFunctionMap map;
map.insert(std::make_pair("Hat", boost::bind( &Mantid::MDAlgorithms::SmoothMD::hatSmooth, instance, _1, _2 )));
map.insert(std::make_pair("Hat", boost::bind( &Mantid::MDAlgorithms::SmoothMD::hatSmooth, instance, _1, _2, _3 )));
return map;
}

Expand Down Expand Up @@ -93,11 +102,13 @@ const std::string SmoothMD::summary() const {
* width.
* @param toSmooth : Workspace to smooth
* @param widthVector : Width vector
* @param weightingWS : Weighting workspace (optional)
* @return Smoothed MDHistoWorkspace
*/
IMDHistoWorkspace_sptr SmoothMD::hatSmooth(IMDHistoWorkspace_const_sptr toSmooth,
const WidthVector &widthVector) {
const WidthVector &widthVector, OptionalIMDHistoWorkspace_const_sptr weightingWS) {

const bool useWeights = weightingWS.is_initialized();
uint64_t nPoints = toSmooth->getNPoints();
Progress progress(this, 0, 1, size_t( double(nPoints) * 1.1 ) );
// Create the output workspace.
Expand All @@ -118,31 +129,53 @@ IMDHistoWorkspace_sptr SmoothMD::hatSmooth(IMDHistoWorkspace_const_sptr toSmooth

do {
// Gets all vertex-touching neighbours
size_t iteratorIndex = iterator->getLinearIndex();

if(useWeights) {

// Check that we could measuer here.
if ( (*weightingWS)->getSignalAt(iteratorIndex) == 0 ) {

outWS->setSignalAt(iteratorIndex, std::numeric_limits<double>::quiet_NaN());

outWS->setErrorSquaredAt(iteratorIndex, std::numeric_limits<double>::quiet_NaN());

continue; // Skip we couldn't measure here.
}
}

std::vector<size_t> neighbourIndexes =
iterator->findNeighbourIndexesByWidth(
widthVector.front()); // TODO we should use the whole width vector
// not just the first element.
const size_t nNeighbours = neighbourIndexes.size();
size_t nNeighbours = neighbourIndexes.size();
double sumSignal = iterator->getSignal();
double sumSqError = iterator->getError();
for (size_t i = 0; i < neighbourIndexes.size(); ++i) {
if(useWeights) {
if ( (*weightingWS)->getSignalAt(neighbourIndexes[i]) == 0 )
{
// Nothing measured here. We cannot use that neighbouring point.
nNeighbours -= 1;
continue;
}
}
sumSignal += toSmooth->getSignalAt(neighbourIndexes[i]);
double error = toSmooth->getErrorAt(neighbourIndexes[i]);
sumSqError += (error * error);

}

// Calculate the mean
outWS->setSignalAt(iterator->getLinearIndex(),
outWS->setSignalAt(iteratorIndex,
sumSignal / double(nNeighbours + 1));
// Calculate the sample variance
outWS->setErrorSquaredAt(iterator->getLinearIndex(),
outWS->setErrorSquaredAt(iteratorIndex,
sumSqError / double(nNeighbours + 1));


progress.report();


} while (iterator->next());
PARALLEL_END_INTERUPT_REGION
}
Expand Down Expand Up @@ -185,6 +218,10 @@ void SmoothMD::init() {
Direction::Input),
docBuffer.str());

declareProperty(new WorkspaceProperty<API::IMDHistoWorkspace>(
"InputNormalizationWorkspace", "", Direction::Input, PropertyMode::Optional),
"Multidimensional weighting workspace. Optional.");

declareProperty(new WorkspaceProperty<API::IMDHistoWorkspace>(
"OutputWorkspace", "", Direction::Output),
"An output smoothed MDHistoWorkspace.");
Expand All @@ -197,6 +234,15 @@ void SmoothMD::exec() {

// Get the input workspace to smooth
IMDHistoWorkspace_sptr toSmooth = this->getProperty("InputWorkspace");

// Get the input weighting workspace
IMDHistoWorkspace_sptr weightingWS = this->getProperty("InputNormalizationWorkspace");
OptionalIMDHistoWorkspace_const_sptr optionalWeightingWS;
if(weightingWS)
{
optionalWeightingWS = weightingWS;
}

// Get the width vector
const std::vector<int> widthVector = this->getProperty("WidthVector");

Expand All @@ -205,7 +251,7 @@ void SmoothMD::exec() {
SmoothFunctionMap functionMap = makeFunctionMap(this);
SmoothFunction smoothFunction = functionMap[smoothFunctionName];
// invoke the smoothing operation
auto smoothed = smoothFunction(toSmooth, widthVector);
auto smoothed = smoothFunction(toSmooth, widthVector, optionalWeightingWS);

setProperty("OutputWorkspace", smoothed);
}
Expand All @@ -218,6 +264,7 @@ std::map<std::string, std::string> SmoothMD::validateInputs() {

std::map<std::string, std::string> product;

// Check the width vector
const std::string widthVectorPropertyName = "WidthVector";
std::vector<int> widthVector = this->getProperty(widthVectorPropertyName);
for (auto it = widthVector.begin(); it != widthVector.end(); ++it) {
Expand All @@ -229,6 +276,43 @@ std::map<std::string, std::string> SmoothMD::validateInputs() {
product.insert(std::make_pair(widthVectorPropertyName, message.str()));
}
}

// Check the dimensionality of the normalization workspace
const std::string normalisationWorkspacePropertyName = "InputNormalizationWorkspace";
IMDHistoWorkspace_sptr toSmoothWs = this->getProperty("InputWorkspace");
IMDHistoWorkspace_sptr normWs = this->getProperty(normalisationWorkspacePropertyName);
if(normWs)
{
const size_t nDimsNorm = normWs->getNumDims();
const size_t nDimsSmooth = toSmoothWs->getNumDims();
if(nDimsNorm != nDimsSmooth) {
std::stringstream message;
message << normalisationWorkspacePropertyName
<< " has a different number of dimensions than InputWorkspace. Shapes of inputs must be the same. Cannot continue smoothing.";
product.insert(std::make_pair(normalisationWorkspacePropertyName, message.str()));
}
else
{
// Loop over dimensions and check nbins.
for(size_t i = 0; i < nDimsNorm; ++i)
{
const size_t nBinsNorm = normWs->getDimension(i)->getNBins();
const size_t nBinsSmooth = toSmoothWs->getDimension(i)->getNBins();
if(nBinsNorm != nBinsSmooth)
{
std::stringstream message;
message << normalisationWorkspacePropertyName
<< ". Number of bins from dimension with index "<< i << " do not match. "
<< nBinsSmooth << " expected. Got " << nBinsNorm
<< ". Shapes of inputs must be the same. Cannot continue smoothing.";
product.insert(std::make_pair(normalisationWorkspacePropertyName, message.str()));
break;
}
}
}
}


return product;
}

Expand Down

0 comments on commit ca69654

Please sign in to comment.