diff --git a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt index 6e3b46216f86a8ed650ecc49f645cf55bb0ccd94..613549cab8e555bf37c2692532f7b19b08a498fe 100644 --- a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt +++ b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt @@ -34,6 +34,7 @@ set ( SRC_FILES src/FitMW.cpp src/FitPowderDiffPeaks.cpp src/FlatBackground.cpp + src/FullprofPolynomial.cpp src/GSLFunctions.cpp src/GSLMatrix.cpp src/GausDecay.cpp @@ -133,6 +134,7 @@ set ( INC_FILES inc/MantidCurveFitting/FitMW.h inc/MantidCurveFitting/FitPowderDiffPeaks.h inc/MantidCurveFitting/FlatBackground.h + inc/MantidCurveFitting/FullprofPolynomial.h inc/MantidCurveFitting/GSLFunctions.h inc/MantidCurveFitting/GSLJacobian.h inc/MantidCurveFitting/GSLMatrix.h @@ -223,6 +225,7 @@ set ( TEST_FILES FitMWTest.h FitPowderDiffPeaksTest.h FlatBackgroundTest.h + FullprofPolynomialTest.h FunctionFactoryConstraintTest.h GSLMatrixTest.h GausDecayTest.h diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/FullprofPolynomial.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/FullprofPolynomial.h new file mode 100644 index 0000000000000000000000000000000000000000..ade3e5d8d8f08008c56f6e44bd02edebb2e5c911 --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/FullprofPolynomial.h @@ -0,0 +1,87 @@ +#ifndef MANTID_CURVEFITTING_FULLPROFPOLYNOMIAL_H_ +#define MANTID_CURVEFITTING_FULLPROFPOLYNOMIAL_H_ + +#include "MantidKernel/System.h" +#include "MantidCurveFitting/BackgroundFunction.h" + +namespace Mantid +{ +namespace CurveFitting +{ + + /** FullprofPolynomial : Polynomial background defined in Fullprof + + Copyright © 2013 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 <http://www.gnu.org/licenses/>. + + File change history is stored at: <https://github.com/mantidproject/mantid> + Code Documentation is available at: <http://doxygen.mantidproject.org> + */ + class DLLExport FullprofPolynomial : public BackgroundFunction + { + public: + FullprofPolynomial(); + virtual ~FullprofPolynomial(); + + /// Overwrite IFunction base class + std::string name()const{return "FullprofPolynomial";} + + virtual const std::string category() const { return "Background";} + + virtual void function1D(double* out, const double* xValues, const size_t nData)const; + + virtual void functionDeriv1D(API::Jacobian* out, const double* xValues, const size_t nData); + + // virtual void functionLocal(std::vector<double> &out, std::vector<double> xValues) const; + + /// Returns the number of attributes associated with the function (polynomial order n) + size_t nAttributes()const{return 1;} + + /// Returns a list of attribute names + std::vector<std::string> getAttributeNames()const; + + /// Return a value of attribute attName + Attribute getAttribute(const std::string& attName)const; + + /// Set a value to attribute attName + void setAttribute(const std::string& attName,const Attribute& ); + + /// Check if attribute attName exists + bool hasAttribute(const std::string& attName)const; + + private: + + /// Polynomial order + int m_n; + + /// Background origin position + double m_bkpos; + + /// Lower x boundary. + double m_StartX; + + /// Upper x boundary + double m_EndX; + }; + + typedef boost::shared_ptr<FullprofPolynomial> FullprofPolynomial_sptr; + + +} // namespace CurveFitting +} // namespace Mantid + +#endif /* MANTID_CURVEFITTING_FULLPROFPOLYNOMIAL_H_ */ diff --git a/Code/Mantid/Framework/CurveFitting/src/FullprofPolynomial.cpp b/Code/Mantid/Framework/CurveFitting/src/FullprofPolynomial.cpp new file mode 100644 index 0000000000000000000000000000000000000000..fc05e841f9cc9bb69f9c42f8b628a40bf5d44fa5 --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/src/FullprofPolynomial.cpp @@ -0,0 +1,186 @@ +#include "MantidCurveFitting/FullprofPolynomial.h" +#include "MantidAPI/FunctionFactory.h" +#include <boost/lexical_cast.hpp> + +namespace Mantid +{ +namespace CurveFitting +{ + + DECLARE_FUNCTION(FullprofPolynomial) + + //---------------------------------------------------------------------------------------------- + /** Constructor + */ + FullprofPolynomial::FullprofPolynomial():m_n(6), m_bkpos(1.) + { + // Declare first 6th order polynomial as default + for(int i=0; i<m_n; ++i) + { + std::string parName = "A" + boost::lexical_cast<std::string>(i); + declareParameter(parName); + } + } + + //---------------------------------------------------------------------------------------------- + /** Destructor + */ + FullprofPolynomial::~FullprofPolynomial() + { + } + + //---------------------------------------------------------------------------------------------- + /** Function to calcualteFullprofPolynomial + */ + void FullprofPolynomial::function1D(double* out, const double* xValues, const size_t nData)const + { + // Generate a vector for all coefficient + std::vector<double> B(m_n, 0.0); + for (int i = 0; i < m_n; ++i) + B[i] = getParameter(i); + + // Calculate + for (size_t i = 0; i < nData; ++i) + { + double tof = xValues[i]; + double x = tof/m_bkpos - 1.; + double pow_x =1.; +#if 1 + // It is the first try as benchmark + double y_b = 0.; + for (int j = 0; j < m_n; ++j) + { + y_b += B[j]*pow_x; + pow_x *= x; + } +#else + // This is the efficient coding + double y_b = B[0]; + for (size_t j = 1; j < m_n; ++j) + { + pow_x *= x; + y_b += B[j] * pow_x; + } +#endif + + out[i] = y_b; + } + + return; + } + + //---------------------------------------------------------------------------------------------- + /** Function to calculate derivative analytically + */ + void FullprofPolynomial::functionDeriv1D(API::Jacobian* out, const double* xValues, const size_t nData) + { + for (size_t i = 0; i < nData; i++) + { + double tof = xValues[i]; + double x = tof/m_bkpos - 1.; + + // Order 0 + double pow_x = 1.0; + out->set(i, 0, pow_x); + // Rest + for (int j = 1; j < m_n; ++j) + { + // It does dirivative to B_j + pow_x *= x; + out->set(i, j, pow_x); + } + } + + return; + } + + //---------------------------------------------------------------------------------------------- + /** Get Attribute names + * @return A list of attribute names (identical toFullprofPolynomial) + */ + std::vector<std::string> FullprofPolynomial::getAttributeNames()const + { + std::vector<std::string> res; + res.push_back("n"); + res.push_back("Bkpos"); + + return res; + } + + //---------------------------------------------------------------------------------------------- + /** Get Attribute + * @param attName :: Attribute name. If it is not "n" exception is thrown. + * @return a value of attribute attName + * (identical toFullprofPolynomial) + */ + API::IFunction::Attribute FullprofPolynomial::getAttribute(const std::string& attName)const + { + if (attName == "n") + { + return Attribute(m_n); + } + else if (attName == "Bkpos") + return Attribute(m_bkpos); + + throw std::invalid_argument("Polynomial: Unknown attribute " + attName); + } + + //---------------------------------------------------------------------------------------------- + /** Set Attribute + * @param attName :: The attribute name. If it is not "n" exception is thrown. + * @param att :: An int attribute containing the new value. The value cannot be negative. + * (identical toFullprofPolynomial) + */ + void FullprofPolynomial::setAttribute(const std::string& attName,const API::IFunction::Attribute& att) + { + if (attName == "n") + { + // set theFullprofPolynomial order + int attint = att.asInt(); + if (attint < 0) + { + throw std::invalid_argument("Polynomial:FullprofPolynomial order cannot be negative."); + } + else if (attint != 6 && attint != 12) + { + throw std::runtime_error("FullprofPolynomial's order must be either 6 or 12. "); + } + else if (attint != m_n) + { + // Only order is (either 6 or 12) and different from current order + clearAllParameters(); + + m_n = attint; + for(int i=0; i<m_n; ++i) + { + std::string parName = "A" + boost::lexical_cast<std::string>(i); + declareParameter(parName); + } + } + } + else if (attName == "Bkpos") + { + // Background original position + m_bkpos = att.asDouble(); + } + + return; + } + + //---------------------------------------------------------------------------------------------- + /** Check if attribute attName exists + */ + bool FullprofPolynomial::hasAttribute(const std::string& attName)const + { + bool has = false; + if (attName == "n") + has = true; + else if (attName == "Bkpos") + has = true; + + return has; + } + + +} // namespace CurveFitting +} // namespace Mantid diff --git a/Code/Mantid/Framework/CurveFitting/test/FullprofPolynomialTest.h b/Code/Mantid/Framework/CurveFitting/test/FullprofPolynomialTest.h new file mode 100644 index 0000000000000000000000000000000000000000..72db146dc0884a1c35d45e079dd03ed32d7cb8d5 --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/test/FullprofPolynomialTest.h @@ -0,0 +1,123 @@ +#ifndef MANTID_CURVEFITTING_FULLPROFPOLYNOMIALTEST_H_ +#define MANTID_CURVEFITTING_FULLPROFPOLYNOMIALTEST_H_ + +#include <cxxtest/TestSuite.h> + +#include "MantidCurveFitting/FullprofPolynomial.h" +#include "MantidAPI/IFunction.h" +#include "MantidAPI/WorkspaceFactory.h" +#include "MantidCurveFitting/Fit.h" +#include "MantidDataObjects/Workspace2D.h" + + +using Mantid::CurveFitting::FullprofPolynomial; + +using namespace Mantid; +using namespace Mantid::API; +using namespace Mantid::CurveFitting; +using namespace Mantid::DataObjects; +using namespace Mantid::Kernel; + +class FullprofPolynomialTest : 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 FullprofPolynomialTest *createSuite() { return new FullprofPolynomialTest(); } + static void destroySuite( FullprofPolynomialTest *suite ) { delete suite; } + + + void testForCategories() + { + FullprofPolynomial forCat; + const std::vector<std::string> categories = forCat.categories(); + TS_ASSERT( categories.size() == 1 ); + TS_ASSERT( categories[0] == "Background" ); + } + + /** Test function on a Fullprof polynomial function + */ + void test_FPPolynomial() + { + // Create a workspace + std::string wsName = "TOFPolybackgroundBackgroundTest"; + int histogramNumber = 1; + int timechannels = 1000; + DataObjects::Workspace2D_sptr ws2D = + boost::dynamic_pointer_cast<DataObjects::Workspace2D>( + API::WorkspaceFactory::Instance().create("Workspace2D", histogramNumber, timechannels, timechannels)); + + AnalysisDataService::Instance().add(wsName, ws2D); + + double tof0 = 8000.; + double dtof = 5.; + + for (int i = 0; i < timechannels; i++) + { + ws2D->dataX(0)[i] = static_cast<double>(i) * dtof + tof0; + } + + // Create a function + IFunction_sptr tofbkgd = boost::dynamic_pointer_cast<IFunction>(boost::make_shared<FullprofPolynomial>()); + TS_ASSERT_THROWS_NOTHING(tofbkgd->setAttributeValue("n", 6)); + TS_ASSERT_THROWS_NOTHING(tofbkgd->setAttributeValue("Bkpos", 10000.)); + TS_ASSERT_THROWS_NOTHING(tofbkgd->setParameter("A0", 0.3)); + TS_ASSERT_THROWS_NOTHING(tofbkgd->setParameter("A1", 1.0)); + TS_ASSERT_THROWS_NOTHING(tofbkgd->setParameter("A2", -0.5)); + TS_ASSERT_THROWS_NOTHING(tofbkgd->setParameter("A3", 0.05)); + TS_ASSERT_THROWS_NOTHING(tofbkgd->setParameter("A4", -0.02)); + + // Calculate function + FunctionDomain1DVector domain(ws2D->readX(0)); + FunctionValues values(domain); + tofbkgd->function(domain, values); + + // Test result + TS_ASSERT_DELTA(values[400], 0.3, 1.0E-10); // Y[10000] = B0 + TS_ASSERT_DELTA(values[0], 0.079568, 1.0E-5); + TS_ASSERT_DELTA(values[605], 0.39730, 1.0E-5); + TS_ASSERT_DELTA(values[999], 0.55583, 1.0E-5); + + // Set the workspace + for (size_t i = 0; i < ws2D->readY(0).size(); ++i) + { + ws2D->dataY(0)[i] = values[i]; + ws2D->dataE(0)[i] = sqrt(fabs(values[i])); + } + + // Make function a little bit off + tofbkgd->setParameter("A0", 0.5); + tofbkgd->setParameter("A3", 0.0); + + // Set up fit + CurveFitting::Fit fitalg; + TS_ASSERT_THROWS_NOTHING(fitalg.initialize()); + TS_ASSERT( fitalg.isInitialized() ); + + fitalg.setProperty("Function", tofbkgd); + fitalg.setPropertyValue("InputWorkspace", wsName); + fitalg.setPropertyValue("WorkspaceIndex", "0"); + + // execute fit + TS_ASSERT_THROWS_NOTHING(TS_ASSERT( fitalg.execute() )); + TS_ASSERT( fitalg.isExecuted() ); + + // test the output from fit is what you expect + double chi2 = fitalg.getProperty("OutputChi2overDoF"); + + TS_ASSERT_DELTA( chi2, 0.0, 0.1); + TS_ASSERT_DELTA( tofbkgd->getParameter("A0"), 0.3, 0.01); + TS_ASSERT_DELTA( tofbkgd->getParameter("A1"), 1.0, 0.0003); + TS_ASSERT_DELTA( tofbkgd->getParameter("A3"), 0.05, 0.01); + + // Clean + AnalysisDataService::Instance().remove(wsName); + + return; + } + + +}; + + +#endif /* MANTID_CURVEFITTING_FULLPROFPOLYNOMIALTEST_H_ */