-
Campbell, Stuart authoredCampbell, Stuart authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
SliceMD.cpp 12.75 KiB
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidKernel/System.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidMDAlgorithms/SliceMD.h"
#include "MantidGeometry/MDGeometry/MDImplicitFunction.h"
#include "MantidKernel/ThreadScheduler.h"
#include "MantidKernel/ThreadPool.h"
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidAPI/FileProperty.h"
#include "MantidGeometry/MDGeometry/MDHistoDimension.h"
#include "MantidKernel/BoundedValidator.h"
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using namespace Mantid::DataObjects;
namespace Mantid {
namespace MDAlgorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(SliceMD)
//----------------------------------------------------------------------------------------------
/** Constructor
*/
SliceMD::SliceMD() {}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
SliceMD::~SliceMD() {}
//----------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void SliceMD::init() {
declareProperty(new WorkspaceProperty<IMDWorkspace>("InputWorkspace", "",
Direction::Input),
"An input MDWorkspace.");
// Properties for specifying the slice to perform.
this->initSlicingProps();
declareProperty(new WorkspaceProperty<Workspace>("OutputWorkspace", "",
Direction::Output),
"Name of the output MDEventWorkspace.");
std::vector<std::string> exts;
exts.push_back(".nxs");
declareProperty(
new FileProperty("OutputFilename", "", FileProperty::OptionalSave, exts),
"Optional: Specify a NeXus file to write if you want the output "
"workspace to be file-backed.");
declareProperty(
new PropertyWithValue<int>("Memory", -1),
"If OutputFilename is specified to use a file back end:\n"
" The amount of memory (in MB) to allocate to the in-memory cache.\n"
" If not specified, a default of 40% of free physical memory is used.");
// setPropertySettings("Memory", new EnabledWhenProperty("OutputFilename",
// IS_NOT_DEFAULT));
declareProperty("TakeMaxRecursionDepthFromInput", true,
"Copy the maximum recursion depth from the input workspace.");
auto mustBePositiveInteger = boost::make_shared<BoundedValidator<int>>();
mustBePositiveInteger->setLower(0);
declareProperty("MaxRecursionDepth", 1000, mustBePositiveInteger,
"Sets the maximum recursion depth to use. Can be used to "
"constrain the workspaces internal structure");
setPropertySettings("MaxRecursionDepth",
new EnabledWhenProperty("TakeMaxRecursionDepthFromInput",
IS_EQUAL_TO, "0"));
setPropertyGroup("OutputFilename", "File Back-End");
setPropertyGroup("Memory", "File Back-End");
}
//----------------------------------------------------------------------------------------------
/** Copy the extra data (not signal, error or coordinates) from one event to
*another
* with different numbers of dimensions
*
* @param srcEvent :: the source event, being copied
* @param newEvent :: the destination event
*/
template <size_t nd, size_t ond>
inline void copyEvent(const MDLeanEvent<nd> &srcEvent,
MDLeanEvent<ond> &newEvent) {
// Nothing extra copy - this is no-op
UNUSED_ARG(srcEvent);
UNUSED_ARG(newEvent);
}
//----------------------------------------------------------------------------------------------
/** Copy the extra data (not signal, error or coordinates) from one event to
*another
* with different numbers of dimensions
*
* @param srcEvent :: the source event, being copied
* @param newEvent :: the destination event
*/
template <size_t nd, size_t ond>
inline void copyEvent(const MDEvent<nd> &srcEvent, MDEvent<ond> &newEvent) {
newEvent.setDetectorId(srcEvent.getDetectorID());
newEvent.setRunIndex(srcEvent.getRunIndex());
}
//----------------------------------------------------------------------------------------------
/** Perform the slice from nd input dimensions to ond output dimensions
*
* @param ws :: input workspace with nd dimensions
* @tparam OMDE :: MDEvent type for the OUTPUT workspace
* @tparam ond :: number of dimensions in the OUTPUT workspace
*/
template <typename MDE, size_t nd, typename OMDE, size_t ond>
void SliceMD::slice(typename MDEventWorkspace<MDE, nd>::sptr ws) {
// Create the ouput workspace
typename MDEventWorkspace<OMDE, ond>::sptr outWS(
new MDEventWorkspace<OMDE, ond>());
for (size_t od = 0; od < m_binDimensions.size(); od++) {
outWS->addDimension(m_binDimensions[od]);
}
outWS->setCoordinateSystem(ws->getSpecialCoordinateSystem());
outWS->initialize();
// Copy settings from the original box controller
BoxController_sptr bc = ws->getBoxController();
// store wrute buffer size for the future
// uint64_t writeBufSize =
// bc->getFileIO()getDiskBuffer().getWriteBufferSize();
// and disable write buffer (if any) for input MD Events for this algorithm
// purposes;
// bc->setCacheParameters(1,0);
BoxController_sptr obc = outWS->getBoxController(); // Use the "number of bins" as the "split into" parameter
for (size_t od = 0; od < m_binDimensions.size(); od++)
obc->setSplitInto(od, m_binDimensions[od]->getNBins());
obc->setSplitThreshold(bc->getSplitThreshold());
bool bTakeDepthFromInputWorkspace =
getProperty("TakeMaxRecursionDepthFromInput");
int tempDepth = getProperty("MaxRecursionDepth");
size_t maxDepth =
bTakeDepthFromInputWorkspace ? bc->getMaxDepth() : size_t(tempDepth);
obc->setMaxDepth(maxDepth);
// size_t outputSize = writeBufSize;
// obc->setCacheParameters(sizeof(OMDE),outputSize);
obc->resetNumBoxes();
// Perform the first box splitting
outWS->splitBox();
size_t lastNumBoxes = obc->getTotalNumMDBoxes();
// --- File back end ? ----------------
std::string filename = getProperty("OutputFilename");
if (!filename.empty()) {
// First save to the NXS file
g_log.notice() << "Running SaveMD to create file back-end" << std::endl;
IAlgorithm_sptr alg = createChildAlgorithm("SaveMD");
alg->setPropertyValue("Filename", filename);
alg->setProperty("InputWorkspace", outWS);
alg->setProperty("MakeFileBacked", true);
alg->executeAsChildAlg();
if (!obc->isFileBacked())
throw std::runtime_error("SliceMD with file-backed output: Can not set "
"up file-backed output workspace ");
auto IOptr = obc->getFileIO();
size_t outBufSize = IOptr->getWriteBufferSize();
// the buffer size for resulting workspace; reasonable size is at least 10
// data chunk sizes (nice to verify)
if (outBufSize < 10 * IOptr->getDataChunk()) {
outBufSize = 10 * IOptr->getDataChunk();
IOptr->setWriteBufferSize(outBufSize);
}
}
// Function defining which events (in the input dimensions) to place in the
// output
MDImplicitFunction *function = this->getImplicitFunctionForChunk(NULL, NULL);
std::vector<API::IMDNode *> boxes;
// Leaf-only; no depth limit; with the implicit function passed to it.
ws->getBox()->getBoxes(boxes, 1000, true, function);
// Sort boxes by file position IF file backed. This reduces seeking time,
// hopefully.
bool fileBackedWS = bc->isFileBacked();
if (fileBackedWS)
API::IMDNode::sortObjByID(boxes);
Progress *prog = new Progress(this, 0.0, 1.0, boxes.size());
// The root of the output workspace
MDBoxBase<OMDE, ond> *outRootBox = outWS->getBox();
// if target workspace has events, we should count them as added
uint64_t totalAdded = outWS->getNEvents();
uint64_t numSinceSplit = 0;
// Go through every box for this chunk.
// PARALLEL_FOR_IF( !bc->isFileBacked() ) for (int i = 0; i < int(boxes.size()); i++) {
MDBox<MDE, nd> *box = dynamic_cast<MDBox<MDE, nd> *>(boxes[i]);
// Perform the binning in this separate method.
if (box) {
// An array to hold the rotated/transformed coordinates
coord_t outCenter[ond];
const std::vector<MDE> &events = box->getConstEvents();
typename std::vector<MDE>::const_iterator it = events.begin();
typename std::vector<MDE>::const_iterator it_end = events.end();
for (; it != it_end; it++) {
// Cache the center of the event (again for speed)
const coord_t *inCenter = it->getCenter();
if (function->isPointContained(inCenter)) {
// Now transform to the output dimensions
m_transformFromOriginal->apply(inCenter, outCenter);
// Create the event
OMDE newEvent(it->getSignal(), it->getErrorSquared(), outCenter);
// Copy extra data, if any
copyEvent(*it, newEvent);
// Add it to the workspace
outRootBox->addEvent(newEvent);
numSinceSplit++;
}
}
box->releaseEvents();
// Ask BC if one needs to split boxes
if (obc->shouldSplitBoxes(totalAdded, numSinceSplit, lastNumBoxes))
// if (numSinceSplit > 20000000 || (i == int(boxes.size()-1)))
{
// This splits up all the boxes according to split thresholds and sizes.
Kernel::ThreadScheduler *ts = new ThreadSchedulerFIFO();
ThreadPool tp(ts);
outWS->splitAllIfNeeded(ts);
tp.joinAll();
// Accumulate stats
totalAdded += numSinceSplit;
numSinceSplit = 0;
lastNumBoxes = obc->getTotalNumMDBoxes();
// Progress reporting
if (!fileBackedWS)
prog->report(i);
}
if (fileBackedWS) {
if (!(i % 10))
prog->report(i);
}
} // is box
} // for each box in the vector
prog->report();
outWS->splitAllIfNeeded(NULL);
// Refresh all cache.
outWS->refreshCache();
g_log.notice() << totalAdded << " " << OMDE::getTypeName()
<< "'s added to the output workspace." << std::endl;
if (outWS->isFileBacked()) {
// Update the file-back-end
g_log.notice() << "Running SaveMD" << std::endl;
IAlgorithm_sptr alg = createChildAlgorithm("SaveMD");
alg->setProperty("UpdateFileBackEnd", true);
alg->setProperty("InputWorkspace", outWS); alg->executeAsChildAlg();
}
// Pass on the display normalization from the input event workspace to the
// output event workspace
IMDEventWorkspace_sptr outEvent =
boost::dynamic_pointer_cast<IMDEventWorkspace>(outWS);
outEvent->setDisplayNormalization(ws->displayNormalization());
outEvent->setDisplayNormalizationHisto(ws->displayNormalizationHisto());
// return the size of the input workspace write buffer to its initial value
// bc->setCacheParameters(sizeof(MDE),writeBufSize);
this->setProperty("OutputWorkspace", outEvent);
delete prog;
}
//----------------------------------------------------------------------------------------------
/// Helper method
template <typename MDE, size_t nd>
void SliceMD::doExec(typename MDEventWorkspace<MDE, nd>::sptr ws) {
if (m_outD == 0)
throw std::runtime_error("No output dimensions specified!");
// Templated method needs to call another templated method depending on the #
// of output dimensions.
if (MDE::getTypeName() == "MDLeanEvent") {
if (m_outD == 1)
this->slice<MDE, nd, MDLeanEvent<1>, 1>(ws);
else if (m_outD == 2)
this->slice<MDE, nd, MDLeanEvent<2>, 2>(ws);
else if (m_outD == 3)
this->slice<MDE, nd, MDLeanEvent<3>, 3>(ws);
else if (m_outD == 4)
this->slice<MDE, nd, MDLeanEvent<4>, 4>(ws);
else
throw std::runtime_error(
"Number of output dimensions > 4. This is not currently handled.");
} else if (MDE::getTypeName() == "MDEvent") {
if (m_outD == 1)
this->slice<MDE, nd, MDEvent<1>, 1>(ws);
else if (m_outD == 2)
this->slice<MDE, nd, MDEvent<2>, 2>(ws);
else if (m_outD == 3)
this->slice<MDE, nd, MDEvent<3>, 3>(ws);
else if (m_outD == 4)
this->slice<MDE, nd, MDEvent<4>, 4>(ws);
else
throw std::runtime_error(
"Number of output dimensions > 4. This is not currently handled.");
} else
throw std::runtime_error("Unexpected MDEvent type '" + MDE::getTypeName() +
"'. This is not currently handled.");
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void SliceMD::exec() {
// Input MDEventWorkspace
m_inWS = getProperty("InputWorkspace");
// Run through the properties to create the transform you need
createTransform();
CALL_MDEVENT_FUNCTION(this->doExec, m_inWS);
}
} // namespace Mantid
} // namespace DataObjects