Skip to content

Commit

Permalink
refs #11393. Cache iterator calculations.
Browse files Browse the repository at this point in the history
Cache the permutation info by width.
  • Loading branch information
OwenArnold committed Mar 20, 2015
1 parent 44c41f5 commit 95407c0
Show file tree
Hide file tree
Showing 3 changed files with 121 additions and 45 deletions.
Expand Up @@ -6,11 +6,14 @@
#include "MantidMDEvents/MDHistoWorkspace.h"
#include "MantidGeometry/MDGeometry/MDImplicitFunction.h"
#include "MantidMDEvents/SkippingPolicy.h"
#include <set>
#include <map>

namespace Mantid {
namespace MDEvents {

// Typdef for a map for mapping width of neighbours (key) to permutations needed in the calcualtion.
typedef std::map<size_t, std::vector<int64_t> > PermutationsMap;

/** An implementation of IMDIterator that iterates through
a MDHistoWorkspace. It treats the bin in the workspace as
a box containing a single "event", at the center of each bin
Expand Down Expand Up @@ -114,6 +117,8 @@ class DLLExport MDHistoWorkspaceIterator : public Mantid::API::IMDIterator {

virtual bool isWithinBounds(size_t index) const;

size_t permutationCacheSize() const;

protected:
/// The MDHistoWorkspace being iterated.
const MDHistoWorkspace *m_ws;
Expand Down Expand Up @@ -153,10 +158,17 @@ class DLLExport MDHistoWorkspaceIterator : public Mantid::API::IMDIterator {

/// Neighbour finding permutations.
mutable std::vector<int64_t> m_permutationsVertexTouching;
/// Neigbour finding permutations for face touching neighbours (3 by 3 width).
mutable std::vector<int64_t> m_permutationsFaceTouching;

/// Neighbour finding permutations map for vertex touching. Keyed via the width (n-pixels) of neighbours required.
mutable PermutationsMap m_permutationsVertexTouchingMap;

/// Skipping policy.
SkippingPolicy_scptr m_skippingPolicy;

/// Create or fetch permutations relating to a given neighbour width.
std::vector<int64_t> createPermutations(const int width) const;
};

} // namespace MDEvents
Expand Down
130 changes: 86 additions & 44 deletions Code/Mantid/Framework/MDEvents/src/MDHistoWorkspaceIterator.cpp
Expand Up @@ -4,6 +4,7 @@
#include "MantidKernel/Utils.h"
#include "MantidGeometry/MDGeometry/IMDDimension.h"
#include <algorithm>
#include <utility>

using namespace Mantid::Kernel;
using namespace Mantid::API;
Expand Down Expand Up @@ -470,15 +471,29 @@ MDHistoWorkspaceIterator::findNeighbourIndexesFaceTouching() const {
iterator.
*/
bool MDHistoWorkspaceIterator::isWithinBounds(size_t index) const {
return index >= m_begin && index < m_max;
return index >= m_begin && index < m_max;
}


std::vector<size_t> MDHistoWorkspaceIterator::findNeighbourIndexesByWidth(const int width) const
{
if( width % 2 == 0 )
{
throw std::invalid_argument("MDHistoWorkspaceIterator::findNeighbourIndexesByWidth, width must always be an even number");
/**
* This is to create the permutations needed to operate find neighbours in the
*vertex-touching schenarios
* Rather than having to fabricate what the possible permutations are each time
*the iterator is moved and the method is called,
* we can cache the results, and re-use them as the only factors are the and the
*dimensionality, the width (n-neighbours).
*
* @return
*/
std::vector<int64_t>
MDHistoWorkspaceIterator::createPermutations(const int width) const {
// look-up
PermutationsMap::iterator it = m_permutationsVertexTouchingMap.find(width);
if (it == m_permutationsVertexTouchingMap.end()) {

if (width % 2 == 0) {
throw std::invalid_argument("MDHistoWorkspaceIterator::"
"findNeighbourIndexesByWidth, width must "
"always be an even number");
}

int64_t offset = 1;
Expand All @@ -487,56 +502,83 @@ std::vector<size_t> MDHistoWorkspaceIterator::findNeighbourIndexesByWidth(const
std::vector<int64_t> permutationsVertexTouching;
permutationsVertexTouching.reserve(integerPower(width, m_nd));

int centreIndex = (width)/2; // Deliberately truncate to get centre index
int centreIndex = (width) / 2; // Deliberately truncate to get centre index

for(int i = 0; i < width; ++i)
{
// for width = 3 : -1, 0, 1
// for width = 5 : -2, -1, 0, 1, 2
permutationsVertexTouching.push_back(centreIndex - i);
for (int i = 0; i < width; ++i) {
// for width = 3 : -1, 0, 1
// for width = 5 : -2, -1, 0, 1, 2
permutationsVertexTouching.push_back(centreIndex - i);
}

// Figure out what possible indexes deltas to generate indexes that are next
// to the current one.
for (size_t j = 1; j < m_nd; ++j) {
offset = offset * static_cast<int64_t>(m_ws->getDimension(j - 1)->getNBins());


size_t nEntries = permutationsVertexTouching.size();
for(size_t k = 1; k <= width/2 ; ++k)
{
for(size_t m = 0; m < nEntries ; m++ )
{
permutationsVertexTouching.push_back( (offset * k) + permutationsVertexTouching[m] );
permutationsVertexTouching.push_back( (offset * k * (-1)) + permutationsVertexTouching[m] );
}
}
offset =
offset * static_cast<int64_t>(m_ws->getDimension(j - 1)->getNBins());

size_t nEntries = permutationsVertexTouching.size();
for (size_t k = 1; k <= width / 2; ++k) {
for (size_t m = 0; m < nEntries; m++) {
permutationsVertexTouching.push_back((offset * k) +
permutationsVertexTouching[m]);
permutationsVertexTouching.push_back((offset * k * (-1)) +
permutationsVertexTouching[m]);
}
}
}

m_permutationsVertexTouchingMap.insert(std::make_pair(width, permutationsVertexTouching));
}

Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker,
m_indexMax, m_index);
// In either case, get the result.
return m_permutationsVertexTouchingMap[width];
}

// Filter out indexes that are are not actually neighbours.
// Accumulate neighbour indexes.
std::vector<size_t> neighbourIndexes;
for (size_t i = 0; i < permutationsVertexTouching.size(); ++i) {
if (permutationsVertexTouching[i] == 0) {
continue;
}
/**
* Find vertex-touching neighbours.
* @param width : Odd number of pixels for all dimensions.
* @return collection of indexes.
*/
std::vector<size_t>
MDHistoWorkspaceIterator::findNeighbourIndexesByWidth(const int width) const {

size_t neighbour_index = m_pos + permutationsVertexTouching[i];
if (neighbour_index < m_ws->getNPoints() &&
Utils::isNeighbourOfSubject(m_nd, neighbour_index, m_index,
m_indexMaker, m_indexMax, width/2)) {
neighbourIndexes.push_back(neighbour_index);
}
// Find existing or create required index permutations.
std::vector<int64_t> permutationsVertexTouching = createPermutations(width);

Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker,
m_indexMax, m_index);

// Filter out indexes that are are not actually neighbours.
// Accumulate neighbour indexes.
std::vector<size_t> neighbourIndexes;
for (size_t i = 0; i < permutationsVertexTouching.size(); ++i) {
if (permutationsVertexTouching[i] == 0) {
continue;
}

size_t neighbour_index = m_pos + permutationsVertexTouching[i];
if (neighbour_index < m_ws->getNPoints() &&
Utils::isNeighbourOfSubject(m_nd, neighbour_index, m_index,
m_indexMaker, m_indexMax, width / 2)) {
neighbourIndexes.push_back(neighbour_index);
}
}

// Remove duplicates
std::sort(neighbourIndexes.begin(), neighbourIndexes.end());
neighbourIndexes.erase(
std::unique(neighbourIndexes.begin(), neighbourIndexes.end()),
neighbourIndexes.end());
return neighbourIndexes;
}

// Remove duplicates
std::sort( neighbourIndexes.begin(), neighbourIndexes.end() );
neighbourIndexes.erase( std::unique( neighbourIndexes.begin(), neighbourIndexes.end() ), neighbourIndexes.end() );
return neighbourIndexes;
/**
*
* @return The size of the permutation cache.
*/
size_t Mantid::MDEvents::MDHistoWorkspaceIterator::permutationCacheSize() const
{
return m_permutationsVertexTouchingMap.size();
}

} // namespace Mantid
Expand Down
22 changes: 22 additions & 0 deletions Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceIteratorTest.h
Expand Up @@ -864,6 +864,28 @@ class MDHistoWorkspaceIteratorTest: public CxxTest::TestSuite
TS_ASSERT(doesContainIndex(neighbourIndexes, 26));
}

void test_cache()
{
const size_t nd = 1;
MDHistoWorkspace_sptr ws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0, nd, 10);
/*
1D MDHistoWorkspace
0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9
*/

MDHistoWorkspaceIterator * it = new MDHistoWorkspaceIterator(ws);
TSM_ASSERT_EQUALS("Empty cache expected", 0, it->permutationCacheSize());
it->findNeighbourIndexesByWidth(3);
TSM_ASSERT_EQUALS("One cache item expected", 1, it->permutationCacheSize());
it->findNeighbourIndexesByWidth(3);
TSM_ASSERT_EQUALS("One cache item expected", 1, it->permutationCacheSize()); // Same item, no change to cache
it->findNeighbourIndexesByWidth(5);
TSM_ASSERT_EQUALS("Two cache entries expected", 2, it->permutationCacheSize());
}


};

class MDHistoWorkspaceIteratorTestPerformance: public CxxTest::TestSuite
Expand Down

0 comments on commit 95407c0

Please sign in to comment.