Newer
Older
#ifndef MANTID_ALGORITHMS_CALCULATEPOLYNOMIALBACKGROUNDTEST_H_
#define MANTID_ALGORITHMS_CALCULATEPOLYNOMIALBACKGROUNDTEST_H_
#include <cxxtest/TestSuite.h>
#include "MantidAlgorithms/CalculatePolynomialBackground.h"
#include "MantidAPI/FrameworkManager.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidHistogramData/Histogram.h"
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
using namespace Mantid;
class CalculatePolynomialBackgroundTest : public CxxTest::TestSuite {
public:
// This pair of boilerplate methods prevent the suite being created statically
// This means the constructor isn't called when running other tests
static CalculatePolynomialBackgroundTest *createSuite() {
return new CalculatePolynomialBackgroundTest();
}
static void destroySuite(CalculatePolynomialBackgroundTest *suite) {
delete suite;
}
CalculatePolynomialBackgroundTest() { API::FrameworkManager::Instance(); }
void test_Init() {
Algorithms::CalculatePolynomialBackground alg;
alg.setRethrows(true);
TS_ASSERT_THROWS_NOTHING(alg.initialize())
TS_ASSERT(alg.isInitialized())
void test_successfulExecutionWithDefaultParameters() {
using namespace WorkspaceCreationHelper;
const auto nHist = 2;
const auto nBin = 2;
auto ws = create2DWorkspaceWhereYIsWorkspaceIndex(nHist, nBin + 1);
auto alg = makeAlgorithm();
TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
TS_ASSERT_THROWS_NOTHING(alg->execute())
TS_ASSERT(alg->isExecuted())
void test_constantBackground() {
using namespace WorkspaceCreationHelper;
const size_t nHist{2};
const size_t nBin{3};
auto ws = create2DWorkspaceWhereYIsWorkspaceIndex(nHist, nBin);
for (size_t histI = 0; histI < nHist; ++histI) {
ws->setCountVariances(histI, nBin, static_cast<double>(histI + 1));
}
auto alg = makeAlgorithm();
TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("Degree", 0))
TS_ASSERT_THROWS_NOTHING(alg->execute())
TS_ASSERT(alg->isExecuted())
API::MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
TS_ASSERT(outWS)
for (size_t histI = 0; histI < nHist; ++histI) {
const auto &ys = ws->y(histI);
const auto &xs = ws->x(histI);
const auto &bkgYs = outWS->y(histI);
const auto &bkgEs = outWS->e(histI);
const auto &bkgXs = outWS->x(histI);
for (size_t binI = 0; binI < nBin; ++binI) {
TS_ASSERT_DELTA(bkgYs[binI], ys[binI], 1e-12)
TS_ASSERT_EQUALS(bkgEs[binI], 0)
TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
}
}
void test_linearBackground() {
using namespace WorkspaceCreationHelper;
const size_t nHist{2};
const size_t nBin{3};
auto ws = create2DWorkspaceWhereYIsWorkspaceIndex(nHist, nBin);
for (size_t histI = 0; histI < nHist; ++histI) {
ws->setCountVariances(histI, nBin, static_cast<double>(histI + 1));
}
auto alg = makeAlgorithm();
TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("Degree", 1))
TS_ASSERT_THROWS_NOTHING(alg->execute())
TS_ASSERT(alg->isExecuted())
API::MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
TS_ASSERT(outWS)
for (size_t histI = 0; histI < nHist; ++histI) {
const auto &ys = ws->y(histI);
const auto &xs = ws->x(histI);
const auto &bkgYs = outWS->y(histI);
const auto &bkgEs = outWS->e(histI);
const auto &bkgXs = outWS->x(histI);
for (size_t binI = 0; binI < nBin; ++binI) {
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
TS_ASSERT_DELTA(bkgYs[binI], ys[binI], 1e-10)
TS_ASSERT_EQUALS(bkgEs[binI], 0)
TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
}
}
}
void test_costFuctionLeastSquares() {
using namespace WorkspaceCreationHelper;
const size_t nHist{2};
const HistogramData::Counts counts{0, 4, 0, 0};
const HistogramData::CountStandardDeviations stdDevs{0, 0.001, 0, 0};
const HistogramData::BinEdges edges{0, 1, 2, 3, 4};
API::MatrixWorkspace_sptr ws(DataObjects::create<DataObjects::Workspace2D>(nHist, HistogramData::Histogram(edges, counts, stdDevs)).release());
auto alg = makeAlgorithm();
TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("Degree", 0))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("CostFunction", "Least squares"))
TS_ASSERT_THROWS_NOTHING(alg->execute())
TS_ASSERT(alg->isExecuted())
API::MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
TS_ASSERT(outWS)
for (size_t histI = 0; histI < nHist; ++histI) {
const auto &xs = ws->x(histI);
const auto &bkgYs = outWS->y(histI);
const auto &bkgEs = outWS->e(histI);
const auto &bkgXs = outWS->x(histI);
for (size_t binI = 0; binI < counts.size(); ++binI) {
// Number 4 in counts is heavily weighted by the small error.
TS_ASSERT_DELTA(bkgYs[binI], 4, 1e-5)
TS_ASSERT_EQUALS(bkgEs[binI], 0)
TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
}
}
}
void test_costFuctionUnweightedLeastSquares() {
using namespace WorkspaceCreationHelper;
const size_t nHist{2};
const HistogramData::Counts counts{0, 4, 0, 0};
const HistogramData::BinEdges edges{0, 1, 2, 3, 4};
API::MatrixWorkspace_sptr ws(DataObjects::create<DataObjects::Workspace2D>(nHist, HistogramData::Histogram(edges, counts)).release());
auto alg = makeAlgorithm();
TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("Degree", 0))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("CostFunction", "Unweighted least squares"))
TS_ASSERT_THROWS_NOTHING(alg->execute())
TS_ASSERT(alg->isExecuted())
API::MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
TS_ASSERT(outWS)
// Unweighted fitting of zeroth order polynomial is equivalent to the mean.
const double result = std::accumulate(counts.cbegin(), counts.cend(), 0.0) / static_cast<double>(counts.size());
for (size_t histI = 0; histI < nHist; ++histI) {
const auto &xs = ws->x(histI);
const auto &bkgYs = outWS->y(histI);
const auto &bkgEs = outWS->e(histI);
const auto &bkgXs = outWS->x(histI);
for (size_t binI = 0; binI < counts.size(); ++binI) {
TS_ASSERT_DELTA(bkgYs[binI], result, 1e-5)
TS_ASSERT_EQUALS(bkgEs[binI], 0)
TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
}
}
}
void test_cubicBackgroundWithNoisyData() {
const double xMin{1000.0};
const double xMax{5000.0};
const double xStep{10.0};
const auto nBins = static_cast<size_t>((xMax - xMin) / xStep);
HistogramData::BinEdges edges(nBins + 1);
{
auto &bins = edges.mutableRawData();
for (size_t i = 0; i < bins.size(); ++i) {
bins[i] = xMin + xStep * static_cast<double>(i);
}
}
HistogramData::Counts counts(nBins);
{
auto &ys = counts.mutableRawData();
for (size_t i = 0; i < ys.size(); ++i) {
// The noise is not random, but a high frequency sinusoidal wave.
const auto x = edges[i] + xStep / 2;
ys[i] = 1000.0 + std::sin(x / 1000.0) + 0.5 * std::sin(x / 40.0);
}
}
HistogramData::Histogram h{edges, counts};
API::MatrixWorkspace_sptr ws =
DataObjects::create<DataObjects::Workspace2D>(1, h);
auto alg = makeAlgorithm();
TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("Degree", 3))
TS_ASSERT_THROWS_NOTHING(alg->execute())
TS_ASSERT(alg->isExecuted())
API::MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
TS_ASSERT(outWS)
TS_ASSERT_EQUALS(outWS->getNumberHistograms(), 1)
const auto &outH = outWS->histogram(0);
TS_ASSERT_EQUALS(outH.size(), h.size())
for (size_t i = 0; i < h.size(); ++i) {
TS_ASSERT_EQUALS(outH.x()[i], h.x()[i])
const auto x = h.x()[i] + xStep / 2;
// Now without the "noise".
const auto y = 1000.0 + std::sin(x / 1000.0);
const auto diff = std::abs(outH.y()[i] - y);
TS_ASSERT_LESS_THAN(diff, 0.08)
TS_ASSERT_EQUALS(outH.e()[i], 0)
}
}
void test_rangesWithGap() {
using namespace WorkspaceCreationHelper;
const size_t nHist{1};
const HistogramData::BinEdges edges{0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5};
const auto nBin = edges.size() - 1;
const HistogramData::Counts counts{1.0, 2.0, 0.0, 0.0, 5.0, 6.0};
const HistogramData::Histogram h{edges, counts};
auto ws = API::MatrixWorkspace_sptr(
DataObjects::create<DataObjects::Workspace2D>(nHist, h));
auto alg = makeAlgorithm();
TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
TS_ASSERT_THROWS_NOTHING(alg->setProperty("Degree", 1))
const std::vector<double> ranges{0.0, 2.5, 4.5, 7.0};
TS_ASSERT_THROWS_NOTHING(alg->setProperty("XRanges", ranges))
TS_ASSERT_THROWS_NOTHING(alg->execute())
TS_ASSERT(alg->isExecuted())
API::MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
TS_ASSERT(outWS)
const auto &xs = ws->x(0);
const auto &bkgYs = outWS->y(0);
const auto &bkgEs = outWS->e(0);
const auto &bkgXs = outWS->x(0);
const std::vector<double> expected{1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
for (size_t binI = 0; binI < nBin; ++binI) {
TS_ASSERT_DELTA(bkgYs[binI], expected[binI], 1e-10)
TS_ASSERT_EQUALS(bkgEs[binI], 0)
TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
}
}
private:
static boost::shared_ptr<Algorithms::CalculatePolynomialBackground>
makeAlgorithm() {
auto a = boost::make_shared<Algorithms::CalculatePolynomialBackground>();
a->initialize();
a->setChild(true);
a->setRethrows(true);
return a;
}
class CalculatePolynomialBackgroundTestPerformance : public CxxTest::TestSuite {
public:
void setUp() override {
Mantid::API::FrameworkManager::Instance();
const size_t nBins = 1000;
const size_t nHisto = 50000;
const std::vector<double> y(nBins, 1.0);
Mantid::HistogramData::Counts counts(y);
std::vector<double> x(nBins + 1);
std::iota(x.begin(), x.end(), 0.0);
const Mantid::HistogramData::BinEdges edges(x);
const Mantid::HistogramData::Histogram h(edges, counts);
m_ws = Mantid::DataObjects::create<Mantid::DataObjects::Workspace2D>(nHisto,
h);
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
// The histograms in m_ws share the same Y and E values. We want unshared
// histograms to performance-test possible cache trashing issues.
for (size_t i = 0; i < m_ws->getNumberHistograms(); ++i) {
const std::vector<double> data(m_ws->y(i).size(), 1.0);
const Mantid::HistogramData::Counts y{data};
const Mantid::HistogramData::CountStandardDeviations e{data};
m_ws->setHistogram(i, edges, y, e);
}
}
void test_zerothDegreePolynomial() {
Mantid::Algorithms::CalculatePolynomialBackground alg;
alg.initialize();
alg.setChild(true);
alg.setProperty("InputWorkspace", m_ws);
alg.setProperty("OutputWorkspace", "__unused_because_child");
alg.setProperty("Degree", 0);
alg.execute();
}
void test_thirdDegreePolynomial() {
Mantid::Algorithms::CalculatePolynomialBackground alg;
alg.initialize();
alg.setChild(true);
alg.setProperty("InputWorkspace", m_ws);
alg.setProperty("OutputWorkspace", "__unused_because_child");
alg.setProperty("Degree", 3);
alg.execute();
}
private:
Mantid::API::MatrixWorkspace_sptr m_ws;
};
#endif /* MANTID_ALGORITHMS_CALCULATEPOLYNOMIALBACKGROUNDTEST_H_ */