diff --git a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt index 809327748560..44b1e2f8c321 100644 --- a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt +++ b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt @@ -67,6 +67,7 @@ set ( SRC_FILES src/NormaliseByPeakArea.cpp src/PRConjugateGradientMinimizer.cpp src/ParDomain.cpp + src/PearsonVIIFamily.cpp src/PlotPeakByLogValue.cpp src/Polynomial.cpp src/ProcessBackground.cpp @@ -177,6 +178,7 @@ set ( INC_FILES inc/MantidCurveFitting/NormaliseByPeakArea.h inc/MantidCurveFitting/PRConjugateGradientMinimizer.h inc/MantidCurveFitting/ParDomain.h + inc/MantidCurveFitting/PearsonVIIFamily.h inc/MantidCurveFitting/PlotPeakByLogValue.h inc/MantidCurveFitting/Polynomial.h inc/MantidCurveFitting/ProcessBackground.h @@ -276,6 +278,7 @@ set ( TEST_FILES NeutronBk2BkExpConvPVoigtTest.h NormaliseByPeakAreaTest.h PRConjugateGradientTest.h + PearsonVIIFamilyTest.h PlotPeakByLogValueTest.h PolynomialTest.h ProcessBackgroundTest.h diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PearsonVIIFamily.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PearsonVIIFamily.h new file mode 100644 index 000000000000..50b802cbf538 --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PearsonVIIFamily.h @@ -0,0 +1,85 @@ +#ifndef MANTID_CURVEFITTING_PEARSONVIIFAMILY_H_ +#define MANTID_CURVEFITTING_PEARSONVIIFAMILY_H_ + +#include "MantidKernel/System.h" +#include "MantidAPI/IFunction1D.h" +#include "MantidAPI/ParamFunction.h" +#include "MantidAPI/IFunction1DAutoDiff.h" + + +namespace Mantid +{ +namespace CurveFitting +{ + + /** PearsonVIIFamily : TODO: DESCRIPTION + + Copyright © 2014 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + File change history is stored at: + Code Documentation is available at: + */ + +namespace Pearsons { + + class DLLExport PearsonVIIHandCoded : public API::IFunction1D, public API::ParamFunction + { + public: + virtual ~PearsonVIIHandCoded() {} + + std::string name() const { return "PearsonVIIHandCoded"; } + + void function1D(double *out, const double *xValues, const size_t nData) const; + void functionDeriv1D(API::Jacobian *out, const double *xValues, const size_t nData); + + protected: + void init(); + }; + + class DLLExport PearsonVIINumDiff : public API::IFunction1D, public API::ParamFunction + { + public: + virtual ~PearsonVIINumDiff() {} + + std::string name() const { return "PearsonVIINumDiff"; } + + void function1D(double *out, const double *xValues, const size_t nData) const; + + protected: + void init(); + }; + + class DLLExport PearsonVIIAutoDiff : virtual public API::IFunction1DAutoDiff, virtual public API::ParamFunction + { + public: + virtual ~PearsonVIIAutoDiff() {} + + std::string name() const { return "PearsonVIIAutoDiff"; } + + void function1DAutoDiff(const API::FunctionDomain1D &domain, std::vector &y, const API::AutoDiffParameterAdapter ¶meters) const; + + protected: + void init(); + + }; + +} // namespace PearsonVIIFamily +} // namespace CurveFitting +} // namespace Mantid + +#endif /* MANTID_CURVEFITTING_PEARSONVIIFAMILY_H_ */ diff --git a/Code/Mantid/Framework/CurveFitting/src/AutoDiffTestAlg.cpp b/Code/Mantid/Framework/CurveFitting/src/AutoDiffTestAlg.cpp index e14caca3d4a1..62e6a4462c33 100644 --- a/Code/Mantid/Framework/CurveFitting/src/AutoDiffTestAlg.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/AutoDiffTestAlg.cpp @@ -3,6 +3,7 @@ #include "MantidCurveFitting/GaussianAutoDiff.h" #include "MantidCurveFitting/GaussianHandCoded.h" #include "MantidCurveFitting/LorentzianFamily.h" +#include "MantidCurveFitting/PearsonVIIFamily.h" #include "MantidAPI/FunctionDomain1D.h" #include "MantidAPI/FunctionValues.h" #include "MantidCurveFitting/Jacobian.h" @@ -60,7 +61,7 @@ const std::string AutoDiffTestAlg::summary() const { return "Test"; } void AutoDiffTestAlg::init() { declareProperty(new WorkspaceProperty("InputWorkspace","",Direction::Input), "Data with 20 gaussian peaks."); - declareProperty("GaussianParameters", "", "Gaussian parameters (h, x0, s)"); + declareProperty("GaussianParameters", "", "Function parameters"); declareProperty(new WorkspaceProperty("OutputWorkspace","",Direction::Output), "Data with 20 gaussian peaks."); declareProperty("DerivativeType","adept","How to calculate derivatives."); @@ -77,14 +78,17 @@ namespace { if ( type == "adept" ) { - return IFunction_sptr(new Lorentzians::LorentzianAutoDiff); - //return IFunction_sptr(new GaussianAutoDiff); + //return IFunction_sptr(new Lorentzians::LorentzianAutoDiff); + return IFunction_sptr(new Pearsons::PearsonVIIAutoDiff); + return IFunction_sptr(new GaussianAutoDiff); } else if ( type == "num") { - return IFunction_sptr(new Lorentzians::LorentzianNumDiff); - //return IFunction_sptr(new GaussianNumDiff); + //return IFunction_sptr(new Lorentzians::LorentzianNumDiff); + return IFunction_sptr(new Pearsons::PearsonVIINumDiff); + return IFunction_sptr(new GaussianNumDiff); } - return IFunction_sptr(new Lorentzians::LorentzianHandCoded); - //return IFunction_sptr(new GaussianHandCoded); + //return IFunction_sptr(new Lorentzians::LorentzianHandCoded); + return IFunction_sptr(new Pearsons::PearsonVIIHandCoded); + return IFunction_sptr(new GaussianHandCoded); } } @@ -101,10 +105,10 @@ void AutoDiffTestAlg::exec() std::string parameterString = getProperty("GaussianParameters"); std::vector paramValues = parameterValues(parameterString); - if(fitData->getNumberHistograms() > 4) { - IFunction_sptr f = getFunction(m_derType); - f->initialize(); + IFunction_sptr f = getFunction(m_derType); + f->initialize(); + if(fitData->getNumberHistograms() > (f->nParams() + 1)) { for(size_t i = 0; i < paramValues.size(); ++i) { f->setParameter(i, paramValues[i] - 0.001 * paramValues[i]); } @@ -115,7 +119,7 @@ void AutoDiffTestAlg::exec() fitAlgorithm->setProperty("CalcErrors", true); fitAlgorithm->setProperty("Function", f); fitAlgorithm->setProperty("InputWorkspace", fitData); - fitAlgorithm->setProperty("WorkspaceIndex", 4); + fitAlgorithm->setProperty("WorkspaceIndex", static_cast(f->nParams() + 1)); fitAlgorithm->execute(); @@ -134,7 +138,7 @@ void AutoDiffTestAlg::exec() FunctionDomain1DVector x(fitData->readX(0)); FunctionValues y(x); - CurveFitting::Jacobian J(x.size(), 3); + CurveFitting::Jacobian J(x.size(), g->nParams()); Kernel::Timer timerF; @@ -166,33 +170,15 @@ void AutoDiffTestAlg::exec() } } - MantidVec &dfdx0 = t->dataY(1); - const MantidVec &dfdx0r = fitData->readY(1); - for(size_t i = 0; i < x.size(); ++i) { - dfdx0[i] = 1.0 - J.get(i, 1) / dfdx0r[i]; - - if(std::isinf(dfdx0[i])) { - dfdx0[i] = 0.0; - } - } - - MantidVec &dfdh = t->dataY(2); - const MantidVec &dfdhr = fitData->readY(2); - for(size_t i = 0; i < x.size(); ++i) { - dfdh[i] = 1.0 - J.get(i, 0) / dfdhr[i]; - - if(std::isinf(dfdh[i])) { - dfdh[i] = 0.0; - } - } - - MantidVec &dfds = t->dataY(3); - const MantidVec &dfdsr = fitData->readY(3); - for(size_t i = 0; i < x.size(); ++i) { - dfds[i] = 1.0 - J.get(i, 2) / dfdsr[i]; + for(size_t j = 1; (j < f->nParams() + 1); ++j) { + MantidVec &partialDeriv = t->dataY(j); + const MantidVec &reference = fitData->readY(j); + for(size_t i = 0; i < partialDeriv.size(); ++i) { + partialDeriv[i] = 1.0 - J.get(i, j - 1) / reference[i]; - if(std::isinf(dfds[i])) { - dfds[i] = 0.0; + if(std::isinf(partialDeriv[i])) { + partialDeriv[i] = 0.0; + } } } diff --git a/Code/Mantid/Framework/CurveFitting/src/GaussianAutoDiff.cpp b/Code/Mantid/Framework/CurveFitting/src/GaussianAutoDiff.cpp index 302be7f523e8..b6c86c297db1 100644 --- a/Code/Mantid/Framework/CurveFitting/src/GaussianAutoDiff.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/GaussianAutoDiff.cpp @@ -30,8 +30,8 @@ void GaussianAutoDiff::function1DAutoDiff(const FunctionDomain1D &domain, std::v } void GaussianAutoDiff::init() { - declareParameter("Height"); declareParameter("PeakCentre"); + declareParameter("Height"); declareParameter("Sigma"); } diff --git a/Code/Mantid/Framework/CurveFitting/src/GaussianHandCoded.cpp b/Code/Mantid/Framework/CurveFitting/src/GaussianHandCoded.cpp index 7500bd9037f7..b7f47e7da129 100644 --- a/Code/Mantid/Framework/CurveFitting/src/GaussianHandCoded.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/GaussianHandCoded.cpp @@ -43,15 +43,15 @@ using namespace API; double term = (xDiff)/sigma; double expTerm = exp(-0.5 * term * term); - out->set(i, 0, expTerm); - out->set(i, 1, (height*expTerm*(xDiff))/(sigma * sigma)); + out->set(i, 0, (height*expTerm*(xDiff))/(sigma * sigma)); + out->set(i, 1, expTerm); out->set(i, 2, (height*expTerm*term * term)/(sigma)); } } void GaussianHandCoded::init() { - declareParameter("Height"); declareParameter("PeakCentre"); + declareParameter("Height"); declareParameter("Sigma"); } diff --git a/Code/Mantid/Framework/CurveFitting/src/GaussianNumDiff.cpp b/Code/Mantid/Framework/CurveFitting/src/GaussianNumDiff.cpp index fb54d99ab07d..71389e849c43 100644 --- a/Code/Mantid/Framework/CurveFitting/src/GaussianNumDiff.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/GaussianNumDiff.cpp @@ -33,8 +33,8 @@ namespace CurveFitting } void GaussianNumDiff::init() { - declareParameter("Height"); declareParameter("PeakCentre"); + declareParameter("Height"); declareParameter("Sigma"); } diff --git a/Code/Mantid/Framework/CurveFitting/src/LorentzianFamily.cpp b/Code/Mantid/Framework/CurveFitting/src/LorentzianFamily.cpp index 755bae9da903..b267800a3c4c 100644 --- a/Code/Mantid/Framework/CurveFitting/src/LorentzianFamily.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/LorentzianFamily.cpp @@ -32,16 +32,16 @@ void LorentzianHandCoded::functionDeriv1D(Jacobian *out, const double *xValues, double num = diff*diff + s * s; double ssquared = s * s; - out->set(i, 0, ssquared / (num)); - out->set(i, 1, 2 * h * ssquared * diff / (num * num)); + out->set(i, 0, 2 * h * ssquared * diff / (num * num)); + out->set(i, 1, ssquared / (num)); out->set(i, 2, 2 * h * s / num * (1 - ssquared / num)); } } void LorentzianHandCoded::init() { - declareParameter("Height"); declareParameter("Centre"); + declareParameter("Height"); declareParameter("Gamma"); } @@ -60,8 +60,8 @@ void LorentzianNumDiff::function1D(double *out, const double *xValues, const siz void LorentzianNumDiff::init() { - declareParameter("Height"); declareParameter("Centre"); + declareParameter("Height"); declareParameter("Gamma"); } @@ -71,7 +71,7 @@ void LorentzianAutoDiff::function1DAutoDiff(const FunctionDomain1D &domain, std: adept::adouble x0 = parameters.getParameter("Centre"); adept::adouble s = parameters.getParameter("Gamma"); - for(size_t i = 0; i < domain.size(); ++i) { + for(size_t i = 0; i < y.size(); ++i) { adept::adouble diff = domain[i] - x0; y[i] = h * s * s / (diff*diff + s * s); } @@ -80,8 +80,8 @@ void LorentzianAutoDiff::function1DAutoDiff(const FunctionDomain1D &domain, std: void LorentzianAutoDiff::init() { - declareParameter("Height"); declareParameter("Centre"); + declareParameter("Height"); declareParameter("Gamma"); } diff --git a/Code/Mantid/Framework/CurveFitting/src/PearsonVIIFamily.cpp b/Code/Mantid/Framework/CurveFitting/src/PearsonVIIFamily.cpp new file mode 100644 index 000000000000..db5a49d174e6 --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/src/PearsonVIIFamily.cpp @@ -0,0 +1,107 @@ +#include "MantidCurveFitting/PearsonVIIFamily.h" + +namespace Mantid +{ +namespace CurveFitting +{ + +namespace Pearsons { + +using namespace API; + +void PearsonVIIHandCoded::function1D(double *out, const double *xValues, const size_t nData) const +{ + double h = getParameter("Height"); + double x0 = getParameter("Centre"); + double s = getParameter("Gamma"); + double m = getParameter("m"); + + for(size_t i = 0; i < nData; ++i) { + double diff = xValues[i] - x0; + out[i] = h * pow(s, 2.0 * m) / pow(diff*diff * (pow(2.0, 1.0/m) - 1.0) + s * s, m); + } +} + +void PearsonVIIHandCoded::functionDeriv1D(Jacobian *out, const double *xValues, const size_t nData) +{ + double h = getParameter("Height"); + double x0 = getParameter("Centre"); + double s = getParameter("Gamma"); + double m = getParameter("m"); + + for(size_t i = 0; i < nData; ++i) { + double diff = xValues[i] - x0; + double ssquared = s * s; + double factor = pow(2.0, (1.0/m)) - 1.0; + double sraised = pow(s, 2.0*m); + double factoredSquaredDiff = factor*diff*diff; + double diffM = pow(ssquared + factoredSquaredDiff, m); + double diffM1 = pow(ssquared + factoredSquaredDiff, m + 1.0); + + out->set(i, 0, 2.0 * (h*m*sraised*diff*factor)/diffM1); + out->set(i, 1, sraised/diffM); + out->set(i, 2, (2.0*h*m) * (pow(s, 2.0*m - 1.0)/diffM - s * sraised/diffM1)); + out->set(i, 3, (2.0*h*sraised*log(s))/diffM + - h*sraised*(log(ssquared + factoredSquaredDiff)/diffM + - ((factor+1.0)*log(2.0)*diff*diff)/(m*diffM1))); + } +} + +void PearsonVIIHandCoded::init() +{ + declareParameter("Centre"); + declareParameter("Height"); + declareParameter("Gamma"); + declareParameter("m"); +} + + +void PearsonVIINumDiff::function1D(double *out, const double *xValues, const size_t nData) const +{ + double h = getParameter("Height"); + double x0 = getParameter("Centre"); + double s = getParameter("Gamma"); + double m = getParameter("m"); + + for(size_t i = 0; i < nData; ++i) { + double diff = xValues[i] - x0; + out[i] = h * pow(s, 2.0 * m) / pow(diff*diff * (pow(2.0, 1.0/m) - 1.0) + s * s, m); + } +} + +void PearsonVIINumDiff::init() +{ + declareParameter("Centre"); + declareParameter("Height"); + declareParameter("Gamma"); + declareParameter("m"); +} + +void PearsonVIIAutoDiff::function1DAutoDiff(const FunctionDomain1D &domain, std::vector &y, const AutoDiffParameterAdapter ¶meters) const +{ + adept::adouble h = parameters.getParameter("Height"); + adept::adouble x0 = parameters.getParameter("Centre"); + adept::adouble s = parameters.getParameter("Gamma"); + adept::adouble m = parameters.getParameter("m"); + + for(size_t i = 0; i < y.size(); ++i) { + adept::adouble diff = domain[i] - x0; + y[i] = h * pow(s, 2.0 * m) / pow(diff*diff * (pow(2.0, 1.0/m) - 1.0) + s * s, m); + } +} + + +void PearsonVIIAutoDiff::init() +{ + declareParameter("Centre"); + declareParameter("Height"); + declareParameter("Gamma"); + declareParameter("m"); +} + +} + + + +} // namespace CurveFitting +} // namespace Mantid diff --git a/Code/Mantid/Framework/CurveFitting/test/PearsonVIIFamilyTest.h b/Code/Mantid/Framework/CurveFitting/test/PearsonVIIFamilyTest.h new file mode 100644 index 000000000000..9f6bbcd5850d --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/test/PearsonVIIFamilyTest.h @@ -0,0 +1,29 @@ +#ifndef MANTID_CURVEFITTING_PEARSONVIIFAMILYTEST_H_ +#define MANTID_CURVEFITTING_PEARSONVIIFAMILYTEST_H_ + +#include + +#include "MantidCurveFitting/PearsonVIIFamily.h" + +using Mantid::CurveFitting::PearsonVIIFamily; +using namespace Mantid::API; + +class PearsonVIIFamilyTest : 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 PearsonVIIFamilyTest *createSuite() { return new PearsonVIIFamilyTest(); } + static void destroySuite( PearsonVIIFamilyTest *suite ) { delete suite; } + + + void test_Something() + { + TSM_ASSERT( "You forgot to write a test!", 0); + } + + +}; + + +#endif /* MANTID_CURVEFITTING_PEARSONVIIFAMILYTEST_H_ */ \ No newline at end of file