Skip to content

Commit

Permalink
Implement lazy caching in SpectraAxis for looking up a spectrum value.
Browse files Browse the repository at this point in the history
Refs #9435
  • Loading branch information
martyngigg committed May 13, 2014
1 parent 94d6f7a commit 1990471
Show file tree
Hide file tree
Showing 7 changed files with 69 additions and 71 deletions.
5 changes: 3 additions & 2 deletions Code/Mantid/Framework/API/inc/MantidAPI/NumericAxis.h
Expand Up @@ -46,6 +46,9 @@ class MatrixWorkspace;
class MANTID_API_DLL NumericAxis: public Axis
{
public:
/// Find the index of the value in the given set of edges
static size_t indexOfValue(const double value, const std::vector<double> & edges);

NumericAxis(const std::size_t& length);
NumericAxis(const std::vector<double>& centres);
virtual ~NumericAxis(){}
Expand Down Expand Up @@ -75,8 +78,6 @@ class MANTID_API_DLL NumericAxis: public Axis
/// Default constructor
NumericAxis();

/// Find the index of the value in the given set of edges
size_t indexOfValue(const double value, const std::vector<double> & edges) const;
/// A vector holding the centre values.
std::vector<double> m_values;

Expand Down
4 changes: 2 additions & 2 deletions Code/Mantid/Framework/API/inc/MantidAPI/SpectraAxis.h
Expand Up @@ -77,8 +77,8 @@ class MANTID_API_DLL SpectraAxis: public Axis
const SpectraAxis& operator=(const SpectraAxis&);
/// A pointer to the workspace holding the axis
const MatrixWorkspace* const m_parentWS;
/// Effective bin width
double m_halfBinWidth;
/// List of edge values for quick searching of values as if this is binned data
mutable std::vector<double> m_edges;
};

} // namespace API
Expand Down
60 changes: 34 additions & 26 deletions Code/Mantid/Framework/API/src/NumericAxis.cpp
Expand Up @@ -12,6 +12,40 @@ namespace Mantid
namespace API
{

//------------------------------------------------------------------------------
// static members
//------------------------------------------------------------------------------

/**
* @param value A value to find in a bin
* @param edges A list of edges that define the bins
* @return The index of the bin holding the value
*/
size_t NumericAxis::indexOfValue(const double value, const std::vector<double> &edges)
{
if(value < edges.front())
{
throw std::out_of_range("NumericAxis::indexOfValue() - Input lower than first value.");
}
auto it = std::lower_bound(edges.begin(), edges.end(), value);
if (it == edges.end())
{
// Out of range
throw std::out_of_range("NumericAxis::indexOfValue() - Input value out of range");
}
// index of closest edge above value is distance of iterator from start
size_t edgeIndex = std::distance(edges.begin(), it);
// index of bin centre is one less since the first boundary offsets the whole range
// need to protect for case where value equals lowest bin boundary as that will return &
// not 1
if(edgeIndex > 0) return edgeIndex - 1;
else return edgeIndex;
}

//------------------------------------------------------------------------------
// public members
//------------------------------------------------------------------------------

/** Constructor
* @param length size of the numeric axis
*/
Expand Down Expand Up @@ -155,31 +189,5 @@ NumericAxis::NumericAxis()
{
}

/**
* @param value A value to find in a bin
* @param edges A list of edges that define the bins
* @return The index of the bin holding the value
*/
size_t NumericAxis::indexOfValue(const double value, const std::vector<double> &edges) const
{
if(value < edges.front())
{
throw std::out_of_range("NumericAxis::indexOfValue() - Input lower than first value.");
}
auto it = std::lower_bound(edges.begin(), edges.end(), value);
if (it == edges.end())
{
// Out of range
throw std::out_of_range("NumericAxis::indexOfValue() - Input value out of range");
}
// index of closest edge above value is distance of iterator from start
size_t edgeIndex = std::distance(edges.begin(), it);
// index of bin centre is one less since the first boundary offsets the whole range
// need to protect for case where value equals lowest bin boundary as that will return &
// not 1
if(edgeIndex > 0) return edgeIndex - 1;
else return edgeIndex;
}

} // namespace API
} // namespace Mantid
41 changes: 13 additions & 28 deletions Code/Mantid/Framework/API/src/SpectraAxis.cpp
Expand Up @@ -3,6 +3,7 @@
//----------------------------------------------------------------------
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/Unit.h"
Expand All @@ -21,10 +22,9 @@ using std::size_t;
* @param parentWorkspace The workspace to which this axis belongs
*/
SpectraAxis::SpectraAxis(const MatrixWorkspace* const parentWorkspace)
: Axis(), m_parentWS(parentWorkspace), m_halfBinWidth(0.5)
: Axis(), m_parentWS(parentWorkspace), m_edges()
{
this->unit() = boost::make_shared<Kernel::Units::Label>("Spectrum", "");
m_halfBinWidth = 0.5*(this->getMax() - this->getMin())/static_cast<double>(m_parentWS->getNumberHistograms()-1);
}

/** Virtual constructor
Expand Down Expand Up @@ -94,36 +94,21 @@ void SpectraAxis::setValue(const std::size_t& index, const double& value)
*/
size_t SpectraAxis::indexOfValue(const double value) const
{
// try and be smart about how we find this. Most of the time
// the index & value don't differ by much. Start at a value
// slightly lower, 2 back, then the actual value so that most cases
// it should find it in a few jumps.
const size_t nspectra = m_parentWS->getNumberHistograms();
const auto & parentWS = *m_parentWS; // avoid constant dereference

size_t guess = (value > 2.0) ? static_cast<size_t>(value - 2.0) : 0;
if(guess >= nspectra) guess = nspectra/2; // start in the middle

size_t index(guess);
do
if(m_edges.empty()) //lazy-instantiation
{
try
m_edges.resize(m_parentWS->getNumberHistograms() + 1);
const size_t npts = m_edges.size() - 1;
for( size_t i = 0; i < npts - 1; ++i )
{
double curNo = static_cast<double>(parentWS.getSpectrum(index)->getSpectrumNo());
if((curNo - m_halfBinWidth <= value) && (value < curNo + m_halfBinWidth)) // creates an effective binning
return index;
m_edges[i+1] = 0.5*(this->getValue(i) + this->getValue(i+1));
}
catch(std::range_error&)
{
continue;
}
++index;
if(index == nspectra) index = 0; //wrap around to search everythin
}
while (index != guess);
// ends
m_edges[0] = this->getValue(0) - (m_edges[1] - this->getValue(0));
m_edges[npts] = this->getValue(npts-1) +
(this->getValue(npts-1) - m_edges[npts-1]);

// Not found if we reached here
throw std::out_of_range("SpectraAxis::indexOfValue() - Value not found on axis");
}
return NumericAxis::indexOfValue(value, m_edges);
}

/** Returns the spectrum number at the position given (Spectra axis only)
Expand Down
13 changes: 11 additions & 2 deletions Code/Mantid/Framework/API/test/SpectraAxisTest.h
Expand Up @@ -114,9 +114,18 @@ class SpectraAxisTest : public CxxTest::TestSuite
{
for (int i=1; i < 6; ++i)
{
TS_ASSERT_EQUALS(i-1, spectraAxis->indexOfValue(static_cast<double>(i - 0.5)) ); //lower boundary considered in this bin
// centre in this bin
TS_ASSERT_EQUALS(i-1, spectraAxis->indexOfValue(static_cast<double>(i)) );
TS_ASSERT_EQUALS(i-1, spectraAxis->indexOfValue(static_cast<double>(i + 0.49)) );
//value on lower boundary in bin below with exception of first boundary where it is above
if(i==1)
{
TS_ASSERT_EQUALS(i-1, spectraAxis->indexOfValue(static_cast<double>(i - 0.5)));
}
else
{
TS_ASSERT_EQUALS(i-2, spectraAxis->indexOfValue(static_cast<double>(i - 0.5)) );
}
TS_ASSERT_EQUALS(i-1, spectraAxis->indexOfValue(static_cast<double>(i + 0.5)) );
}
}

Expand Down
8 changes: 0 additions & 8 deletions Code/Mantid/Framework/DataHandling/test/SaveSPETest.h
Expand Up @@ -177,14 +177,6 @@ class SaveSPETest : public CxxTest::TestSuite
MatrixWorkspace_sptr setUpWorkspace(const std::string & input, MatrixWorkspace_sptr inputWS)
{
inputWS->getAxis(0)->unit() = Mantid::Kernel::UnitFactory::Instance().create("DeltaE");

// the following is largely about associating detectors with the workspace
for (int j = 0; j < NHIST; ++j)
{
// Just set the spectrum number to match the index
inputWS->getAxis(1)->setValue(j, j+1);
}

AnalysisDataService::Instance().add(input,inputWS);

// Load the instrument data
Expand Down
9 changes: 6 additions & 3 deletions Code/Mantid/Framework/DataObjects/src/Workspace2D.cpp
Expand Up @@ -53,9 +53,6 @@ namespace Mantid
{
m_noVectors = NVectors;
data.resize(m_noVectors);
m_axes.resize(2);
m_axes[0] = new API::RefAxis(XLength, this);
m_axes[1] = new API::SpectraAxis(this);

MantidVecPtr t1,t2;
t1.access().resize(XLength); //this call initializes array to zero
Expand All @@ -74,6 +71,12 @@ namespace Mantid
spec->setSpectrumNo(specid_t(i+1));
spec->setDetectorID(detid_t(i+1));
}

// Add axes that reference the data
m_axes.resize(2);
m_axes[0] = new API::RefAxis(XLength, this);
m_axes[1] = new API::SpectraAxis(this);

}


Expand Down

0 comments on commit 1990471

Please sign in to comment.