Skip to content

Commit

Permalink
refs #12694. Refactor to method and code doc.
Browse files Browse the repository at this point in the history
  • Loading branch information
@OwenArnold
OwenArnold committed Sep 8, 2015
1 parent 194ee8d commit 14011b1
Show file tree
Hide file tree
Showing 2 changed files with 138 additions and 43 deletions.
9 changes: 7 additions & 2 deletions Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/ReplicateMD.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,9 @@ namespace Mantid {
namespace API {
class IMDHistoWorkspace;
}
namespace DataObjects {
class MDHistoWorkspace;
}
namespace MDAlgorithms {

/** ReplicateMD : TODO: DESCRIPTION
Expand Down Expand Up @@ -49,9 +52,11 @@ class MANTID_MDALGORITHMS_DLL ReplicateMD : public API::Algorithm {
virtual std::map<std::string, std::string> validateInputs();

private:
boost::shared_ptr<const Mantid::API::IMDHistoWorkspace>
transposeMD(boost::shared_ptr<Mantid::API::IMDHistoWorkspace>& toTranspose,
boost::shared_ptr<const Mantid::DataObjects::MDHistoWorkspace>
transposeMD(boost::shared_ptr<Mantid::DataObjects::MDHistoWorkspace>& toTranspose,
const std::vector<int> &axes);
boost::shared_ptr<Mantid::DataObjects::MDHistoWorkspace> getDataWorkspace() const;
boost::shared_ptr<Mantid::DataObjects::MDHistoWorkspace> getShapeWorkspace() const;
void init();
void exec();
};
Expand Down
172 changes: 131 additions & 41 deletions Code/Mantid/Framework/MDAlgorithms/src/ReplicateMD.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
#include "MantidMDAlgorithms/ReplicateMD.h"
#include "MantidKernel/Utils.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/Progress.h"
#include "MantidAPI/IMDHistoWorkspace.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidDataObjects/MDHistoWorkspaceIterator.h"
#include "MantidGeometry/MDGeometry/IMDDimension.h"
#include <utility>
Expand Down Expand Up @@ -38,6 +39,15 @@ IMDDimension_const_sptr findMatchingDimension(const IMDHistoWorkspace &toSearch,
return foundDim;
}

/**
* Find the index of the dimension to be replicated. This is any dimension
* in the shape workspace that does not appear in the data workspace (there
* should be only 1)
* or any integrated dimension in the data workspace.
* @param shapeWS : workspace with complete target dimension set
* @param dataWS : source of the missing/integrated dimension
* @return Index in the SHAPE workspace
*/
size_t findReplicationDimension(const IMDHistoWorkspace &shapeWS,
const IMDHistoWorkspace &dataWS) {
size_t index = -1;
Expand All @@ -53,8 +63,12 @@ size_t findReplicationDimension(const IMDHistoWorkspace &shapeWS,
return index;
}



/**
* Determine the axis to use for a transpose operation.
* @param shapeWS : Workspace to tranpose into the coordinate system of
* @param dataWS : Workspace to transpose
* @return Axis required for the transpose command.
*/
std::vector<int> findAxes(const IMDHistoWorkspace &shapeWS,
const IMDHistoWorkspace &dataWS) {

Expand All @@ -68,6 +82,44 @@ std::vector<int> findAxes(const IMDHistoWorkspace &shapeWS,
}
return axes;
}

/**
* Convert a linear index in the shape to a linear index in the DataItem
*
* This is a flattening procedure.
*
* 1) Convert linear index in shape to an n-dimensional set of indexes
* 2) Remove the index corresponding to the missing dimension (we're going to replicate over this)
* 3) Using the n-dimensional indexes, create a linear index in the data
*
* @param nDimsShape : Number of dimensions in the shape
* @param shapeReplicationIndex : Index to replicate over
* @param indexMaxShape : Number of bins for each dim in the shape
* @param indexMakerShape : Index maker for the shape
* @param sourceIndex : Linear index in the shape
* @param indexMakerData : Index maker for the data
* @param nDimsData : Number of dimensions in the data
* @return Linear index in the data
*/
size_t linearIndexToLinearIndex(const size_t &nDimsShape,
const size_t &shapeReplicationIndex,
const std::vector<size_t> &indexMaxShape,
const std::vector<size_t> &indexMakerShape,
const size_t &sourceIndex,
std::vector<size_t> &indexMakerData,
const size_t &nDimsData) {
std::vector<size_t> vecShapeIndexes(nDimsShape);
Utils::NestedForLoop::GetIndicesFromLinearIndex(
nDimsShape, sourceIndex, &indexMakerShape[0], &indexMaxShape[0],
&vecShapeIndexes[0]);

vecShapeIndexes.erase(vecShapeIndexes.begin() + shapeReplicationIndex);

const size_t targetIndex = Utils::NestedForLoop::GetLinearIndex(
nDimsData, &vecShapeIndexes[0], &indexMakerData[0]);

return targetIndex;
}
}

// Register the algorithm into the AlgorithmFactory
Expand Down Expand Up @@ -100,18 +152,29 @@ const std::string ReplicateMD::summary() const {
"replicating along an additional axis";
}

IMDHistoWorkspace_const_sptr
ReplicateMD::transposeMD(IMDHistoWorkspace_sptr &toTranspose,
/**
* Transpose the input data workspace according to the axis provided.
* @param toTranspose Workspace to transpose
* @param axes : target axes indexes
* @return : Transposed workspace.
*/
MDHistoWorkspace_const_sptr
ReplicateMD::transposeMD(MDHistoWorkspace_sptr &toTranspose,
const std::vector<int> &axes) {

auto transposeMD = this->createChildAlgorithm("TransposeMD");
auto transposeMD = this->createChildAlgorithm("TransposeMD", 0.0, 0.5);
transposeMD->setProperty("InputWorkspace", toTranspose);
transposeMD->setProperty("Axes", axes);
transposeMD->execute();
IMDHistoWorkspace_sptr outputWS = transposeMD->getProperty("OutputWorkspace");
return outputWS;
return boost::dynamic_pointer_cast<const MDHistoWorkspace>(outputWS);
}

/**
* Overriden validate inputs
*
* @return map of propery name to problem description for any issues
*/
std::map<std::string, std::string> ReplicateMD::validateInputs() {
std::map<std::string, std::string> errorMap;
IMDHistoWorkspace_sptr shapeWS = this->getProperty("ShapeWorkspace");
Expand Down Expand Up @@ -174,22 +237,52 @@ void ReplicateMD::init() {
"An output workspace with replicated data.");
}


//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
/**
* fetches and down casts
* @return ShapeWorkspace property
*/
void ReplicateMD::exec() {
MDHistoWorkspace_sptr ReplicateMD::getShapeWorkspace() const {
MDHistoWorkspace_sptr shapeWS;
{
IMDHistoWorkspace_sptr temp = this->getProperty("ShapeWorkspace");
shapeWS = boost::dynamic_pointer_cast<MDHistoWorkspace>(temp);
}
IMDHistoWorkspace_sptr dataWS = this->getProperty("DataWorkspace");

return shapeWS;
}

/**
* fetches and down casts
* @return DataWorkspace property
*/
MDHistoWorkspace_sptr ReplicateMD::getDataWorkspace() const {
MDHistoWorkspace_sptr dataWS;
{
IMDHistoWorkspace_sptr temp = this->getProperty("DataWorkspace");
dataWS = boost::dynamic_pointer_cast<MDHistoWorkspace>(temp);
}

return dataWS;
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void ReplicateMD::exec() {

MDHistoWorkspace_sptr shapeWS = getShapeWorkspace();
MDHistoWorkspace_sptr dataWS = getDataWorkspace();

const size_t nDimsShape = shapeWS->getNumDims();
size_t nDimsData = dataWS->getNumDims();

IMDHistoWorkspace_const_sptr transposedDataWS = dataWS;
Progress progress(this, 0.0, 1.0, size_t(shapeWS->getNPoints()));

/*
We transpose the input so that we have the same order of dimensions in data
as in the shape. This is important. If they are not in the same order, then
the linear index -> linear index calculation below will not work correctly.
*/
MDHistoWorkspace_const_sptr transposedDataWS = dataWS;
if (dataWS->getNumDims() == shapeWS->getNumDims()) {
auto axes = findAxes(*shapeWS, *dataWS);
transposedDataWS = transposeMD(dataWS, axes);
Expand All @@ -199,51 +292,48 @@ void ReplicateMD::exec() {
// Set up index maker for shape
std::vector<size_t> indexMakerShape(nDimsShape);
std::vector<size_t> indexMaxShape(nDimsShape);
for(size_t i = 0; i < nDimsShape; ++i){
indexMaxShape[i] = shapeWS->getDimension(i)->getNBins();
for (size_t i = 0; i < nDimsShape; ++i) {
indexMaxShape[i] = shapeWS->getDimension(i)->getNBins();
}
Utils::NestedForLoop::SetUpIndexMaker(nDimsShape, &indexMakerShape[0], &indexMaxShape[0]);
Utils::NestedForLoop::SetUpIndexMaker(nDimsShape, &indexMakerShape[0],
&indexMaxShape[0]);

// Set up index maker for data.
std::vector<size_t> indexMakerData(nDimsData);
std::vector<size_t> indexMaxData(nDimsData);
for(size_t i = 0; i < nDimsData; ++i){
indexMaxData[i] = transposedDataWS->getDimension(i)->getNBins();
for (size_t i = 0; i < nDimsData; ++i) {
indexMaxData[i] = transposedDataWS->getDimension(i)->getNBins();
}
Utils::NestedForLoop::SetUpIndexMaker(nDimsData, &indexMakerData[0], &indexMaxData[0]);
Utils::NestedForLoop::SetUpIndexMaker(nDimsData, &indexMakerData[0],
&indexMaxData[0]);

size_t shapeReplicationIndex = findReplicationDimension(*shapeWS, *transposedDataWS);
// Dimension index into the shape workspace where we will be replicating
const size_t shapeReplicationIndex =
findReplicationDimension(*shapeWS, *transposedDataWS);

// Create the output workspace from the shape.
MDHistoWorkspace_sptr outputWS = shapeWS->clone();
auto outIt = std::unique_ptr<MDHistoWorkspaceIterator>(
dynamic_cast<MDHistoWorkspaceIterator *>(outputWS->createIterator()));

auto dataIt = std::unique_ptr<MDHistoWorkspaceIterator>(
dynamic_cast<MDHistoWorkspaceIterator *>(
transposedDataWS->createIterator()));



do {

const auto sourceIndex = outIt->getLinearIndex();

std::vector<size_t> vecShapeIndexes(nDimsShape);
Utils::NestedForLoop::GetIndicesFromLinearIndex(nDimsShape, sourceIndex, &indexMakerShape[0], &indexMaxShape[0], &vecShapeIndexes[0]);

vecShapeIndexes.erase(vecShapeIndexes.begin() + shapeReplicationIndex);

const size_t targetIndex = Utils::NestedForLoop::GetLinearIndex(nDimsData, &vecShapeIndexes[0],
&indexMakerData[0]);

const size_t targetIndex = linearIndexToLinearIndex(
nDimsShape, shapeReplicationIndex, indexMaxShape, indexMakerShape,
sourceIndex, indexMakerData, nDimsData);

// const auto targetIndex = indexInData(sourceIndex, shapeWS);
dataIt->jumpTo(targetIndex);
outputWS->setSignalAt(sourceIndex, dataIt->getSignal());
outputWS->setSignalAt(sourceIndex,
transposedDataWS->getSignalAt(targetIndex));
outputWS->setErrorSquaredAt(sourceIndex,
dataIt->getError() * dataIt->getError());
outputWS->setNumEventsAt(sourceIndex, dataIt->getNumEvents());
outputWS->setMDMaskAt(sourceIndex, dataIt->getIsMasked());
transposedDataWS->getErrorAt(targetIndex) *
transposedDataWS->getErrorAt(targetIndex));
outputWS->setNumEventsAt(sourceIndex,
transposedDataWS->getNumEventsAt(targetIndex));
outputWS->setMDMaskAt(sourceIndex,
transposedDataWS->getIsMaskedAt(targetIndex));
progress.report();

} while (outIt->next());

Expand Down

0 comments on commit 14011b1

Please sign in to comment.