Skip to content
Snippets Groups Projects
Commit f61bf2a3 authored by Zhou, Wenduo's avatar Zhou, Wenduo
Browse files

Merge pull request #13356 from mantidproject/13335_SplineSmoothing_Alg_Taking_Overtime

Fixes SplineSmoothing bugs and executes the algorithm faster.
It works well with the new feature.  Thus it is closed. 
parents 8624fee2 97897638
No related branches found
No related tags found
No related merge requests found
......@@ -61,7 +61,7 @@ private:
void exec();
/// smooth a single spectrum of the workspace
void smoothSpectrum(int index);
void smoothSpectrum(int index, int maxBreaks);
/// calculate derivatives for a single spectrum
void calculateSpectrumDerivatives(int index, int order);
......@@ -81,7 +81,7 @@ private:
/// choose points to define a spline and smooth the data
void selectSmoothingPoints(API::MatrixWorkspace_const_sptr inputWorkspace,
size_t row);
size_t row, int maxBreaks);
/// calculate the spline based on the smoothing points chosen
void calculateSmoothing(API::MatrixWorkspace_const_sptr inputWorkspace,
......
......@@ -69,6 +69,12 @@ void SplineSmoothing::init() {
errorSizeValidator->setLower(0.0);
declareProperty("Error", 0.05, errorSizeValidator,
"The amount of error we wish to tolerate in smoothing");
auto numOfBreaks = boost::make_shared<BoundedValidator<int>>();
numOfBreaks->setLower(0);
declareProperty("MaxNumberOfBreaks", M_START_SMOOTH_POINTS, numOfBreaks,
"To set the positions of the break-points, default 10 "
"equally spaced real values in interval 0.0 - 1.0");
}
//----------------------------------------------------------------------------------------------
......@@ -80,6 +86,9 @@ void SplineSmoothing::exec() {
int histNo = static_cast<int>(m_inputWorkspace->getNumberHistograms());
int order = static_cast<int>(getProperty("DerivOrder"));
// retrieving number of breaks
int maxBreaks = static_cast<int>(getProperty("MaxNumberOfBreaks"));
m_inputWorkspacePointData = convertBinnedData(m_inputWorkspace);
m_outputWorkspace = setupOutputWorkspace(m_inputWorkspacePointData, histNo);
......@@ -90,7 +99,7 @@ void SplineSmoothing::exec() {
Progress pgress(this, 0.0, 1.0, histNo);
for (int i = 0; i < histNo; ++i) {
smoothSpectrum(i);
smoothSpectrum(i, maxBreaks);
calculateSpectrumDerivatives(i, order);
pgress.report();
}
......@@ -105,13 +114,14 @@ void SplineSmoothing::exec() {
/** Smooth a single spectrum of the input workspace
*
* @param index :: index of the spectrum to smooth
* @param maxBreaks :: the number to set the positions of the break-points
*/
void SplineSmoothing::smoothSpectrum(int index) {
void SplineSmoothing::smoothSpectrum(int index, int maxBreaks) {
m_cspline = boost::make_shared<BSpline>();
m_cspline->setAttributeValue("Uniform", false);
// choose some smoothing points from input workspace
selectSmoothingPoints(m_inputWorkspacePointData, index);
selectSmoothingPoints(m_inputWorkspacePointData, index, maxBreaks);
performAdditionalFitting(m_inputWorkspacePointData, index);
// compare the data set against our spline
......@@ -204,10 +214,9 @@ SplineSmoothing::convertBinnedData(MatrixWorkspace_sptr workspace) {
* @param outputWorkspace :: The output workspace
* @param row :: The row of spectra to use
*/
void
SplineSmoothing::calculateSmoothing(MatrixWorkspace_const_sptr inputWorkspace,
MatrixWorkspace_sptr outputWorkspace,
size_t row) const {
void SplineSmoothing::calculateSmoothing(
MatrixWorkspace_const_sptr inputWorkspace,
MatrixWorkspace_sptr outputWorkspace, size_t row) const {
// define the spline's parameters
const auto &xIn = inputWorkspace->readX(row);
size_t nData = xIn.size();
......@@ -300,20 +309,26 @@ void SplineSmoothing::addSmoothingPoints(const std::set<int> &points,
*
* @param inputWorkspace :: The input workspace containing noisy data
* @param row :: The row of spectra to use
* @param maxBreaks :: the number to set the positions of the break-points
*/
void SplineSmoothing::selectSmoothingPoints(
MatrixWorkspace_const_sptr inputWorkspace, size_t row) {
MatrixWorkspace_const_sptr inputWorkspace, size_t row, int maxBreaks) {
std::set<int> smoothPts;
const auto &xs = inputWorkspace->readX(row);
const auto &ys = inputWorkspace->readY(row);
int xSize = static_cast<int>(xs.size());
// if retrienved value is default zero
if (maxBreaks == 0) {
setProperty("MaxNumberOfBreaks", xs.size());
maxBreaks = getProperty("MaxNumberOfBreaks");
}
// number of points to start with
int numSmoothPts(M_START_SMOOTH_POINTS);
int numSmoothPts(maxBreaks);
// evenly space initial points over data set
int delta = xSize / numSmoothPts;
for (int i = 0; i < xSize; i += delta) {
smoothPts.insert(i);
}
......@@ -322,7 +337,7 @@ void SplineSmoothing::selectSmoothingPoints(
bool resmooth(true);
while (resmooth) {
// if we're using all points then we can't do anything more.
if (smoothPts.size() >= xs.size() - 1)
if (smoothPts.size() > (unsigned)(maxBreaks + 2))
break;
addSmoothingPoints(smoothPts, xs.data(), ys.data());
......
......@@ -8,94 +8,92 @@
using Mantid::CurveFitting::SplineSmoothing;
class SplineSmoothingTest : public CxxTest::TestSuite
{
class SplineSmoothingTest : 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 SplineSmoothingTest *createSuite() { return new SplineSmoothingTest(); }
static void destroySuite( SplineSmoothingTest *suite ) { delete suite; }
//Functor for generating data
struct SplineFunc
{
double operator()(double x, int)
{
return std::sin(x);
}
static SplineSmoothingTest *createSuite() {
return new SplineSmoothingTest();
}
static void destroySuite(SplineSmoothingTest *suite) { delete suite; }
// Functor for generating data
struct SplineFunc {
double operator()(double x, int) { return std::sin(x); }
};
void test_Init()
{
void test_Init() {
SplineSmoothing alg;
TS_ASSERT_THROWS_NOTHING( alg.initialize() )
TS_ASSERT( alg.isInitialized() )
TS_ASSERT_THROWS_NOTHING(alg.initialize())
TS_ASSERT(alg.isInitialized())
}
void testExec()
{
void testExec() {
using namespace Mantid::API;
//number of derivatives and spectra
const int order (2), spectra(2);
// number of derivatives and spectra
const int order(2), spectra(2);
//create a binned workspace
MatrixWorkspace_sptr inputWorkspace = WorkspaceCreationHelper::Create2DWorkspaceFromFunction(SplineFunc(), spectra, 0, 5, 0.02, false);
// create a binned workspace
MatrixWorkspace_sptr inputWorkspace =
WorkspaceCreationHelper::Create2DWorkspaceFromFunction(
SplineFunc(), spectra, 0, 5, 0.02, false);
//setup algorithm
// setup algorithm
SplineSmoothing alg;
runAlgorithm(alg, order, inputWorkspace);
checkOutput(alg);
}
void testExecHistogramData()
{
void testExecHistogramData() {
using namespace Mantid::API;
//number of derivatives and spectra
const int order (2), spectra(1);
// number of derivatives and spectra
const int order(2), spectra(1);
//create a binned workspace
MatrixWorkspace_sptr inputWorkspace = WorkspaceCreationHelper::Create2DWorkspaceFromFunction(SplineFunc(), spectra, 0, 5, 0.02, true);
// create a binned workspace
MatrixWorkspace_sptr inputWorkspace =
WorkspaceCreationHelper::Create2DWorkspaceFromFunction(
SplineFunc(), spectra, 0, 5, 0.02, true);
SplineSmoothing alg;
runAlgorithm(alg, order, inputWorkspace);
checkOutput(alg);
}
void testExecMultipleHistograms()
{
void testExecMultipleHistograms() {
using namespace Mantid::API;
//number of derivatives and spectra
const int order (2), spectra(3);
// number of derivatives and spectra
const int order(2), spectra(3);
//create a binned workspace
MatrixWorkspace_sptr inputWorkspace = WorkspaceCreationHelper::Create2DWorkspaceFromFunction(SplineFunc(), spectra, 0, 5, 0.02, true);
// create a binned workspace
MatrixWorkspace_sptr inputWorkspace =
WorkspaceCreationHelper::Create2DWorkspaceFromFunction(
SplineFunc(), spectra, 0, 5, 0.02, true);
SplineSmoothing alg;
runAlgorithm(alg, order, inputWorkspace);
checkOutput(alg);
}
void checkOutput(const SplineSmoothing& alg) const
{
void checkOutput(const SplineSmoothing &alg) const {
using namespace Mantid::API;
MatrixWorkspace_const_sptr ows = alg.getProperty("OutputWorkspace");
WorkspaceGroup_const_sptr derivs = alg.getProperty("OutputWorkspaceDeriv");
for (size_t i = 0; i < ows->getNumberHistograms(); ++i)
{
MatrixWorkspace_const_sptr derivsWs = boost::dynamic_pointer_cast<const MatrixWorkspace>(derivs->getItem(i));
const auto & xs = ows->readX(i);
const auto & ys = ows->readY(i);
const auto & d1 = derivsWs->readY(0);
const auto & d2 = derivsWs->readY(1);
//check output for consistency
for(size_t j = 0; j < ys.size(); ++j)
{
for (size_t i = 0; i < ows->getNumberHistograms(); ++i) {
MatrixWorkspace_const_sptr derivsWs =
boost::dynamic_pointer_cast<const MatrixWorkspace>(
derivs->getItem(i));
const auto &xs = ows->readX(i);
const auto &ys = ows->readY(i);
const auto &d1 = derivsWs->readY(0);
const auto &d2 = derivsWs->readY(1);
// check output for consistency
for (size_t j = 0; j < ys.size(); ++j) {
TS_ASSERT_DELTA(ys[j], std::sin(xs[j]), 1e-4);
TS_ASSERT_DELTA(d1[j], std::cos(xs[j]), 1e-1);
TS_ASSERT_DELTA(d2[j], -std::sin(xs[j]), 1e-1);
......@@ -103,23 +101,22 @@ public:
}
}
void runAlgorithm(SplineSmoothing& alg, int order, const Mantid::API::MatrixWorkspace_sptr& iws) const
{
void runAlgorithm(SplineSmoothing &alg, int order,
const Mantid::API::MatrixWorkspace_sptr &iws) const {
alg.initialize();
alg.setChild(true);
alg.setPropertyValue("OutputWorkspace", "Anon");
alg.setPropertyValue("OutputWorkspaceDeriv", "AnonDerivs");
TS_ASSERT_THROWS_NOTHING( alg.setProperty("Error", 0.05));
TS_ASSERT_THROWS_NOTHING( alg.setProperty("DerivOrder", order));
TS_ASSERT_THROWS_NOTHING(alg.setProperty("Error", 0.05));
TS_ASSERT_THROWS_NOTHING(alg.setProperty("DerivOrder", order));
TS_ASSERT_THROWS_NOTHING(alg.setProperty("MaxNumberOfBreaks", 50));
alg.setProperty("InputWorkspace", iws);
TS_ASSERT_THROWS_NOTHING( alg.execute() );
TS_ASSERT( alg.isExecuted() );
TS_ASSERT_THROWS_NOTHING(alg.execute());
TS_ASSERT(alg.isExecuted());
}
};
#endif /* MANTID_ALGORITHMS_SPLINETEST_H_ */
......@@ -10,14 +10,30 @@ Description
-----------
The algorithm performs a smoothing of the input data using a cubic
spline. The algorithm takes a 2D workspace and generates a spline for
each of the spectra to approximate a fit of the data.
spline. The algorithm takes a 2D workspace and generates a spline
for each of the spectra to approximate a fit of the data.
Optionally, this algorithm can also calculate the first and second
derivatives of each of the interpolated points as a side product.
Setting the DerivOrder property to zero will force the algorithm to
calculate no derivatives.
The algorithm allows user to include number of breaking points which
will enable the algorithm to execute functions with large number of
spikes or noise. With the control of number of breaking points, user
will be able to execute :ref:`SplineSmoothing <algm-SplineSmoothing-v1>`
algorithm in small period of time. The lower breaking points number is
set the faster algorithm will execute the function but it will not
produce high quality smoothing function. By inserting high number of
breaking points, a detailed smoothing function can also be produced.
Users are also able to successfully apply the algorithm to a workspace
with multiple spectrums, which would generate an output of workspace
with multiple spectrums of :ref:`SplineSmoothing <algm-SplineSmoothing-v1>`
algorithm applied to it. `BSpline <http://www.mantidproject.org/BSpline>`_
can be used to help you understand break-points in further detail.
For Histogram Workspaces
########################
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment