Commit ff9e808a authored by Federico Montesino Pouzols's avatar Federico Montesino Pouzols
Browse files

Merge pull request #15692 from mantidproject/15652_Reduce_box_splitting_in_SliceMD

Improvements to SliceMD algorithm
parents 275c97cf 91c6506b
......@@ -207,10 +207,16 @@ TMDE(std::vector<coord_t> MDEventWorkspace)::estimateResolution() const {
if (numMD[i] > 0)
realDepth = i;
auto splitTop = m_BoxController->getSplitTopInto();
std::vector<coord_t> out;
for (size_t d = 0; d < nd; d++) {
size_t finestSplit = 1;
for (size_t i = 0; i < realDepth; i++)
size_t i = 0;
if (splitTop) {
finestSplit *= splitTop.get()[d];
i = 1;
}
for (; i < realDepth; i++)
finestSplit *= m_BoxController->getSplitInto(d);
Geometry::IMDDimension_const_sptr dim = this->getDimension(d);
// Calculate the bin size at the smallest split amount
......
......@@ -242,8 +242,8 @@ TMDE(size_t MDGridBox)::computeSizesFromSplit() {
splitCumul[d] = tot;
tot *= split[d];
// Length of the side of a box in this dimension
m_SubBoxSize[d] =
double(this->extents[d].getSize()) / static_cast<double>(split[d]);
m_SubBoxSize[d] = static_cast<double>(this->extents[d].getSize()) /
static_cast<double>(split[d]);
// Accumulate the squared diagonal length.
diagSum += m_SubBoxSize[d] * m_SubBoxSize[d];
}
......@@ -476,7 +476,8 @@ TMDE(void MDGridBox)::getBoxes(std::vector<API::IMDNode *> &outBoxes,
Kernel::Utils::NestedForLoop::SetUp(nd, vertexIndex, 0);
// To get indexes in the array of vertexes
size_t vertexIndexMaker[nd];
Kernel::Utils::NestedForLoop::SetUpIndexMaker(nd, vertexIndexMaker, vertices_max);
Kernel::Utils::NestedForLoop::SetUpIndexMaker(nd, vertexIndexMaker,
vertices_max);
// To get indexes in the array of BOXES
size_t boxIndexMaker[nd];
Kernel::Utils::NestedForLoop::SetUpIndexMaker(nd, boxIndexMaker, split);
......@@ -521,8 +522,8 @@ TMDE(void MDGridBox)::getBoxes(std::vector<API::IMDNode *> &outBoxes,
if (i & mask)
vertIndex[d] = 1;
}
size_t linIndex =
Kernel::Utils::NestedForLoop::GetLinearIndex(nd, vertIndex, vertexIndexMaker);
size_t linIndex = Kernel::Utils::NestedForLoop::GetLinearIndex(
nd, vertIndex, vertexIndexMaker);
vertexNeighborsOffsets[i] = linIndex;
}
......@@ -533,8 +534,8 @@ TMDE(void MDGridBox)::getBoxes(std::vector<API::IMDNode *> &outBoxes,
bool allDone = false;
while (!allDone) {
// Find the linear index into the BOXES array.
size_t boxLinearIndex =
Kernel::Utils::NestedForLoop::GetLinearIndex(nd, boxIndex, boxIndexMaker);
size_t boxLinearIndex = Kernel::Utils::NestedForLoop::GetLinearIndex(
nd, boxIndex, boxIndexMaker);
API::IMDNode *box = m_Children[boxLinearIndex];
// std::cout << "Box at " << Strings::join(boxIndex, boxIndex+nd,
......@@ -544,8 +545,8 @@ TMDE(void MDGridBox)::getBoxes(std::vector<API::IMDNode *> &outBoxes,
// Find the linear index of the upper left vertex of the box.
// (note that we're using the VERTEX index maker to find the linear index
// in that LARGER array)
size_t vertLinearIndex =
Kernel::Utils::NestedForLoop::GetLinearIndex(nd, boxIndex, vertexIndexMaker);
size_t vertLinearIndex = Kernel::Utils::NestedForLoop::GetLinearIndex(
nd, boxIndex, vertexIndexMaker);
// OK, now its time to see if the box is touching or contained or out of
// it.
......@@ -856,8 +857,8 @@ TMDE(void MDGridBox)::centerpointBin(MDBin<MDE, nd> &bin,
}
// Increment the counter(s) in the nested for loops.
allDone = Kernel::Utils::NestedForLoop::Increment(nd, counters, counters_max,
counters_min);
allDone = Kernel::Utils::NestedForLoop::Increment(
nd, counters, counters_max, counters_min);
}
}
......@@ -1169,8 +1170,8 @@ TMDE(void MDGridBox)::integrateSphere(API::CoordTransform &radiusTransform,
}
if (!badIndex) {
// Convert to linear index
size_t linearIndex =
Kernel::Utils::NestedForLoop::GetLinearIndex(nd, boxIndex, indexMaker);
size_t linearIndex = Kernel::Utils::NestedForLoop::GetLinearIndex(
nd, boxIndex, indexMaker);
// So we have one more vertex touching this box that is contained in
// the integration volume. Whew!
verticesContained[linearIndex]++;
......@@ -1181,7 +1182,8 @@ TMDE(void MDGridBox)::integrateSphere(API::CoordTransform &radiusTransform,
}
// Increment the counter(s) in the nested for loops.
allDone = Kernel::Utils::NestedForLoop::Increment(nd, vertexIndex, vertices_max);
allDone =
Kernel::Utils::NestedForLoop::Increment(nd, vertexIndex, vertices_max);
}
// OK, we've done all the vertices. Now we go through and check each box.
......@@ -1387,8 +1389,8 @@ TMDE(void MDGridBox)::integrateCylinder(
}
if (!badIndex) {
// Convert to linear index
size_t linearIndex =
Kernel::Utils::NestedForLoop::GetLinearIndex(nd, boxIndex, indexMaker);
size_t linearIndex = Kernel::Utils::NestedForLoop::GetLinearIndex(
nd, boxIndex, indexMaker);
// So we have one more vertex touching this box that is contained in
// the integration volume. Whew!
verticesContained[linearIndex]++;
......@@ -1399,7 +1401,8 @@ TMDE(void MDGridBox)::integrateCylinder(
}
// Increment the counter(s) in the nested for loops.
allDone = Kernel::Utils::NestedForLoop::Increment(nd, vertexIndex, vertices_max);
allDone =
Kernel::Utils::NestedForLoop::Increment(nd, vertexIndex, vertices_max);
}
// OK, we've done all the vertices. Now we go through and check each box.
......@@ -1640,6 +1643,12 @@ TMDE(void MDGridBox)::buildAndAddEventUnsafe(const signal_t Signal,
* */
TMDE(inline size_t MDGridBox)::addEvent(const MDE &event) {
size_t cindex = calculateChildIndex(event);
// We can erroneously get cindex == numBoxes for events which fall on the
// upper boundary of the last child box, so add these events to the last box
if (cindex == numBoxes)
cindex = numBoxes - 1;
if (cindex < numBoxes)
return m_Children[cindex]->addEvent(event);
else
......@@ -1665,6 +1674,12 @@ TMDE(inline size_t MDGridBox)::addEvent(const MDE &event) {
* */
TMDE(inline size_t MDGridBox)::addEventUnsafe(const MDE &event) {
size_t cindex = calculateChildIndex(event);
// We can erroneously get cindex == numBoxes for events which fall on the
// upper boundary of the last child box, so add these events to the last box
if (cindex == numBoxes)
cindex = numBoxes - 1;
if (cindex < numBoxes)
return m_Children[cindex]->addEventUnsafe(event);
else
......@@ -1728,7 +1743,7 @@ TMDE(size_t MDGridBox)::calculateChildIndex(const MDE &event) const {
for (size_t d = 0; d < nd; d++) {
// Accumulate the index
auto offset = event.getCenter(d) - this->extents[d].getMin();
cindex += int(offset / (m_SubBoxSize[d])) * splitCumul[d];
cindex += static_cast<int>(offset / (m_SubBoxSize[d])) * splitCumul[d];
}
return cindex;
}
......
......@@ -439,6 +439,28 @@ public:
TS_ASSERT_DELTA(binSizes[1], 1.0, 1e-6);
}
//-------------------------------------------------------------------------------------
void test_estimateResolution_with_top_level_splitting() {
MDEventWorkspace2Lean::sptr b =
MDEventsTestHelper::makeMDEW<2>(10, 0.0, 10.0);
std::vector<coord_t> binSizes;
// First, before any splitting
binSizes = b->estimateResolution();
TS_ASSERT_EQUALS(binSizes.size(), 2);
TS_ASSERT_DELTA(binSizes[0], 10.0, 1e-6);
TS_ASSERT_DELTA(binSizes[1], 10.0, 1e-6);
auto bc = b->getBoxController();
bc->setSplitTopInto(0, 5);
// Resolution is smaller after splitting
b->splitBox();
binSizes = b->estimateResolution();
TS_ASSERT_EQUALS(binSizes.size(), 2);
TS_ASSERT_DELTA(binSizes[0], 2.0, 1e-6);
TS_ASSERT_DELTA(binSizes[1], 10.0, 1e-6);
}
////-------------------------------------------------------------------------------------
void
......
......@@ -135,9 +135,11 @@ void SliceMD::slice(typename MDEventWorkspace<MDE, nd>::sptr ws) {
// 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());
// Use the "number of bins" as the "split into" parameter for the top level
for (size_t od = 0; od < m_binDimensions.size(); od++) {
obc->setSplitTopInto(od, m_binDimensions[od]->getNBins());
obc->setSplitInto(od, bc->getSplitInto(od));
}
obc->setSplitThreshold(bc->getSplitThreshold());
bool bTakeDepthFromInputWorkspace =
......@@ -228,9 +230,8 @@ void SliceMD::slice(typename MDEventWorkspace<MDE, nd>::sptr ws) {
// Copy extra data, if any
copyEvent(*it, newEvent);
// Add it to the workspace
outRootBox->addEvent(newEvent);
numSinceSplit++;
if (outRootBox->addEvent(newEvent))
numSinceSplit++;
}
}
box->releaseEvents();
......@@ -265,8 +266,10 @@ void SliceMD::slice(typename MDEventWorkspace<MDE, nd>::sptr ws) {
// Refresh all cache.
outWS->refreshCache();
// Account for events that were added after the last split
totalAdded += numSinceSplit;
g_log.notice() << totalAdded << " " << OMDE::getTypeName()
<< "'s added to the output workspace." << std::endl;
<< "s added to the output workspace." << std::endl;
if (outWS->isFileBacked()) {
// Update the file-back-end
......
......@@ -394,12 +394,9 @@ void MantidQwtIMDWorkspaceData::choosePlotAxis()
bool atLeastOneDimNotIntegrated = false;
for(size_t i = 0; i < mdew->getNumDims(); ++i)
{
if( mdew->getDimension(i)->getNBins() == controller->getSplitInto(i))
if(!mdew->getDimension(i)->getIsIntegrated())
{
if(!mdew->getDimension(i)->getIsIntegrated())
{
atLeastOneDimNotIntegrated = true;
}
atLeastOneDimNotIntegrated = true;
}
}
regularBinnedMDWorkspace = atLeastOneDimNotIntegrated;
......
......@@ -38,13 +38,18 @@ Splitting Parameters
####################
The **OutputBins** parameter is interpreted as the "SplitInto" parameter
for each dimension. For instance, if you want the output workspace to
for the top split-level in each dimension. For instance, if you want the output workspace to
split in 2x2x2, you would specify OutputBins="2,2,2".
For 1D slices, it may make sense to specify a SplitInto parameter of 1
in every other dimension - that way, boxes will only be split along the
1D direction.
To force the box structure to match that defined in OutputBins, the
MaxRecursionDepth property can be set to 1. If this is not the case then
boxes will split further if sufficient events fall in the same box. Further
splitting uses the value of "SplitInto" from the InputWorkspace.
Slicing a MDHistoWorkspace
##########################
......
......@@ -70,6 +70,8 @@ MD Algorithms (VATES CLI)
``.vtu`` files. These file types can be loaded into a standalone version
of ParaView.
- PlotMD now plots points at bin centres for MDEventWorkspaces as well as MDHistoWorkspaces.
- SliceMD now reports the correct number of events in the output workspace.
- The size of densely populated, multidimensional MDEventWorkspace slices produced by SliceMD has been greatly reduced by using more sensible box splitting parameters.
Performance
-----------
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment