diff --git a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt index 93f89bd5197a..6d16289ad805 100644 --- a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt +++ b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt @@ -281,6 +281,7 @@ set ( TEST_FILES MuonFInteractionTest.h NeutronBk2BkExpConvPVoigtTest.h NormaliseByPeakAreaTest.h + PeakParametersNumericTest.h PRConjugateGradientTest.h PlotPeakByLogValueTest.h PolynomialTest.h diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PeakParametersNumeric.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PeakParametersNumeric.h index ff5ad56332db..1a6819f18087 100644 --- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PeakParametersNumeric.h +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PeakParametersNumeric.h @@ -14,7 +14,7 @@ class ChebfunBase; /** Implements IPeakFunction's getters and setters for the peak centre, height and -fwhm. +fwhm. The parameters are calculated numerically. Copyright © 2007-8 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge National Laboratory & European Spallation Source @@ -57,6 +57,11 @@ class DLLExport PeakParametersNumeric : public API::IPeakFunction { bool explicitlySet = true); using API::IPeakFunction::setParameter; + /// Get boundaries for an interval within which the peak has + /// significant values. The interval must contain the maximum + /// point and points below the half-maximum on both side + /// from the centre. Ideally the interval should exclude + /// points with values below 1% of the maximum. virtual std::pair getExtent() const = 0; protected: @@ -66,25 +71,35 @@ class DLLExport PeakParametersNumeric : public API::IPeakFunction { std::vector &p, std::vector &a) const; + /// Enumerates possible effects of parameters on the width enum WidthParamType {Linear, Square, Inverse}; void defineHeightParameter(const std::string& parName); void defineCentreParameter(const std::string& parName); void defineWidthParameter(const std::string& parName, WidthParamType wType); - mutable bool m_invalidateCache; - mutable double m_centre; - mutable double m_width; - mutable double m_height; - - mutable double m_start; - mutable double m_end; - mutable boost::shared_ptr m_base; + // setter helpers + /// Index of parameter proportional to the height size_t m_heightIndex; + /// Index of parameter which shifts the centre size_t m_centreIndex; + /// Indices of parameters which affect the width std::vector m_widthIndices; + /// Types of the width parameter effects std::vector m_widthParTypes; + + // cached values + + mutable bool m_invalidateCache; ///< dirty flag + mutable double m_centre; ///< peak centre (the maximum point) + mutable double m_width; ///< FWHM + mutable double m_height; ///< the maximum value + + mutable double m_start; + mutable double m_end; + mutable boost::shared_ptr m_base; ///< function approximator + }; } // namespace CurveFitting diff --git a/Code/Mantid/Framework/CurveFitting/src/PeakParametersNumeric.cpp b/Code/Mantid/Framework/CurveFitting/src/PeakParametersNumeric.cpp index 91747a4a2aec..fbdaae0d728c 100644 --- a/Code/Mantid/Framework/CurveFitting/src/PeakParametersNumeric.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/PeakParametersNumeric.cpp @@ -15,32 +15,43 @@ namespace CurveFitting { using namespace Kernel; using namespace API; +/// Specify a name of a parameter that scales in sync with the peak height. +/// Multiplying this parameter by a factor changes the height by the same +/// factor. +/// @param parName :: A parameter name void PeakParametersNumeric::defineHeightParameter(const std::string &parName) { m_heightIndex = parameterIndex(parName); } +/// Specify a name of a parameter that shifts in sync with the peak centre. +/// Adding a value to this parameter moves the centre by the same amount. +/// @param parName :: A parameter name void PeakParametersNumeric::defineCentreParameter(const std::string &parName) { m_centreIndex = parameterIndex(parName); } +/// Add a name of a parameter that affects the peak width. There can be more +/// than one such parameter. +/// @param parName :: A parameter name +/// @param wType :: The kind of dependency between the width and this parameter. +/// If the width needs to be scaled by factor f +/// - Linear parameter is scaled by the same factor f, +/// - Square parameter is scaled by f^2 +/// - Inverse parameter is scaled by 1/f. void PeakParametersNumeric::defineWidthParameter(const std::string &parName, WidthParamType wType) { m_widthIndices.push_back(parameterIndex(parName)); m_widthParTypes.push_back(wType); } -/** - * Get approximate height of the peak: function value at X0. - */ +/// Calculate the peak height as the extreme value. double PeakParametersNumeric::height() const { updateCache(); return m_height; } -/** - * Set new height of the peak. This method does this approximately. - * @param h :: New value for the height. - */ +/// Set new height of the peak. +/// @param h :: New value for the height. void PeakParametersNumeric::setHeight(const double h) { double h0 = height(); if (h0 == 0.0) { @@ -58,11 +69,14 @@ void PeakParametersNumeric::setHeight(const double h) { setParameter(m_heightIndex, parValue); } +/// Calculate the peak centre as the extreme point. double PeakParametersNumeric::centre() const { updateCache(); return m_centre; } +/// Set new centre position. +/// @param c :: New centre value. void PeakParametersNumeric::setCentre(const double c) { double c0 = centre(); double dc = c - c0; @@ -70,18 +84,16 @@ void PeakParametersNumeric::setCentre(const double c) { setParameter(m_centreIndex, x0); } -/** - * Get approximate peak width. - */ +/// Get the peak width as the distance between the two points +/// of half-maximum on both sides of the centre. double PeakParametersNumeric::fwhm() const { updateCache(); return m_width; } -/** - * Set new peak width approximately. - * @param w :: New value for the width. - */ +/// Set new peak width by scaling parameters specified by calls +/// to defineWidthParameter(...). +/// @param w :: New value for the width. void PeakParametersNumeric::setFwhm(const double w) { double wOld = fwhm(); double factor = w / wOld; @@ -103,64 +115,84 @@ void PeakParametersNumeric::setFwhm(const double w) { } } -/** - * Calculate function value for a single argument. - */ +/// Calculate function value for a single argument. +/// @param x :: Function argument. double PeakParametersNumeric::operator()(double x) const { double y = 0.0; function1D(&y, &x, 1); return y; } +/// Override the base class method to set the dirty flag when any +/// of the parameters changes. void PeakParametersNumeric::setParameter(size_t i, const double &value, bool explicitlySet) { IPeakFunction::setParameter(i, value, explicitlySet); m_invalidateCache = true; } +/// Make function approximator ChebfunBase to approximate the peak +/// on an interval. +/// @param start :: Start of the interval. +/// @param end :: End of the interval. +/// @param p :: Peak values at pointes defined by ChebfunBase. +/// @param a :: Chebyshev expansion coefficients. boost::shared_ptr PeakParametersNumeric::makeBase(double start, double end, std::vector &p, std::vector &a) const { double tolerance = 1e-15; ChebfunBase_sptr base; + // Run bestFit with decreasing tolerance until approximation succeeds + // Approximation of high accuracy can fail in cases of too sharp peaks + // or too large extents. while (tolerance < 1.0) { base = ChebfunBase::bestFit(start, end, *this, p, a, 0.0, tolerance, 100); if (base) return base; tolerance *= 100; } + // If all failed create an approximation with whatever accuracy + // we can get. base = boost::make_shared(8, start, end); p = base->fit(*this); a = base->calcA(p); return base; } +/// Calculate the centre, peak and width if dirty flag has been raised. void PeakParametersNumeric::updateCache() const { if (!m_invalidateCache) return; m_invalidateCache = false; + // Get an interval for the approximation const auto interval = getExtent(); double start = interval.first; double end = interval.second; + // Containers for approximation's values std::vector a, p, ad; bool baseBuilt = false; for (size_t iter = 0; iter < 2; ++iter) { baseBuilt = true; - m_base = makeBase(start, end, p, a); + m_base = makeBase(start, end, p, a); m_base->derivative(a, ad); + // Find the root(s) of the derivative which must give peak centre auto roots = m_base->roots(ad); + // If approximation is bad there can be 0 or more than 1 roots if (roots.empty()) { + // set the centre in the middle of the interval m_centre = (start + end) / 2; } else if (roots.size() == 1) { + // this has to be the correct value m_centre = roots[0]; } else { + // if approximation ascillates find the root with the highest value double maxVal = 0.0; size_t iMax = 0; for (size_t i = 0; i < roots.size(); ++i) { @@ -173,8 +205,12 @@ void PeakParametersNumeric::updateCache() const { m_centre = roots[iMax]; } + // height is peak's value at the centre m_height = (*this)(m_centre); + // At this point we can check if getExtent() gave us a good interval. + // Check the peak values at the ends of the interval are below a certain + // fraction of the height. double h = fabs(m_height) / 8; double dStart = m_centre - start; while (fabs((*this)(start)) > h) { @@ -187,14 +223,19 @@ void PeakParametersNumeric::updateCache() const { baseBuilt = false; } + // If the interval turned out to be too small baseBuilt is false + // and we make another iteration if (baseBuilt) break; } + // The fastest way of shifting the approximation down the y-axis by height/2 a[0] -= m_height / 2; + // Now the roots should give the points of half-maximum auto roots = m_base->roots(a); + // If the approximation isn't perfect consider different possibilities if (roots.empty()) { m_width = (end - start) / 2; } else if (roots.size() == 1) { diff --git a/Code/Mantid/Framework/CurveFitting/test/ChebfunBaseTest.h b/Code/Mantid/Framework/CurveFitting/test/ChebfunBaseTest.h index 060167fb0d20..eafa66f878c2 100644 --- a/Code/Mantid/Framework/CurveFitting/test/ChebfunBaseTest.h +++ b/Code/Mantid/Framework/CurveFitting/test/ChebfunBaseTest.h @@ -1,5 +1,5 @@ -#ifndef CHEBYSHEVTEST_H_ -#define CHEBYSHEVTEST_H_ +#ifndef CHEBFUNBASETEST_H_ +#define CHEBFUNBASETEST_H_ #include @@ -249,4 +249,4 @@ class ChebfunBaseTest : public CxxTest::TestSuite { } }; -#endif /*CHEBYSHEVTEST_H_*/ +#endif /*CHEBFUNBASETEST_H_*/