diff --git a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt index a39e4c433ab8..d6ebca85d114 100644 --- a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt +++ b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt @@ -12,6 +12,7 @@ set ( SRC_FILES src/BivariateNormal.cpp src/Bk2BkExpConvPV.cpp src/BoundaryConstraint.cpp + src/CalculateChiSquared.cpp src/CalculateGammaBackground.cpp src/CalculateMSVesuvio.cpp src/Chebyshev.cpp @@ -53,6 +54,7 @@ set ( SRC_FILES src/Gaussian.cpp src/GaussianComptonProfile.cpp src/GramCharlierComptonProfile.cpp + src/IFittingAlgorithm.cpp src/IkedaCarpenterPV.cpp src/LatticeDomainCreator.cpp src/LatticeFunction.cpp @@ -126,6 +128,7 @@ set ( INC_FILES inc/MantidCurveFitting/BivariateNormal.h inc/MantidCurveFitting/Bk2BkExpConvPV.h inc/MantidCurveFitting/BoundaryConstraint.h + inc/MantidCurveFitting/CalculateChiSquared.h inc/MantidCurveFitting/CalculateGammaBackground.h inc/MantidCurveFitting/CalculateMSVesuvio.h inc/MantidCurveFitting/Chebyshev.h @@ -170,6 +173,7 @@ set ( INC_FILES inc/MantidCurveFitting/Gaussian.h inc/MantidCurveFitting/GaussianComptonProfile.h inc/MantidCurveFitting/GramCharlierComptonProfile.h + inc/MantidCurveFitting/IFittingAlgorithm.h inc/MantidCurveFitting/IkedaCarpenterPV.h inc/MantidCurveFitting/Jacobian.h inc/MantidCurveFitting/LatticeDomainCreator.h @@ -240,6 +244,7 @@ set ( TEST_FILES BivariateNormalTest.h Bk2BkExpConvPVTest.h BoundaryConstraintTest.h + CalculateChiSquaredTest.h CalculateGammaBackgroundTest.h CalculateMSVesuvioTest.h ChebyshevTest.h diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CalculateChiSquared.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CalculateChiSquared.h new file mode 100644 index 000000000000..43694114c164 --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CalculateChiSquared.h @@ -0,0 +1,52 @@ +#ifndef MANTID_CURVEFITTING_CALCULATECHISQUARED_H_ +#define MANTID_CURVEFITTING_CALCULATECHISQUARED_H_ + +#include "MantidKernel/System.h" +#include "MantidCurveFitting/IFittingAlgorithm.h" + +namespace Mantid { +namespace CurveFitting { + +/** + + Calculate chi squared for a function and a data set in a workspace. + + Copyright © 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge + National Laboratory & European Spallation Source + + 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: +*/ +class DLLExport CalculateChiSquared : public IFittingAlgorithm { +public: + CalculateChiSquared(); + virtual ~CalculateChiSquared(); + + virtual const std::string name() const; + virtual int version() const; + virtual const std::string summary() const; + +private: + void initConcrete(); + void execConcrete(); +}; + +} // namespace CurveFitting +} // namespace Mantid + +#endif /* MANTID_CURVEFITTING_CALCULATECHISQUARED_H_ */ diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/Fit.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/Fit.h index 754872b5cdd5..7bd44423f203 100644 --- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/Fit.h +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/Fit.h @@ -4,10 +4,11 @@ //---------------------------------------------------------------------- // Includes //---------------------------------------------------------------------- -#include "MantidAPI/Algorithm.h" -#include "MantidAPI/IFunction.h" -#include "MantidAPI/Workspace.h" -#include "MantidAPI/IDomainCreator.h" +//#include "MantidAPI/Algorithm.h" +//#include "MantidAPI/IFunction.h" +//#include "MantidAPI/Workspace.h" +//#include "MantidAPI/IDomainCreator.h" +#include "MantidCurveFitting/IFittingAlgorithm.h" namespace Mantid { @@ -93,10 +94,10 @@ along with this program. If not, see . File change history is stored at: Code Documentation is available at: */ -class DLLExport Fit : public API::Algorithm { +class DLLExport Fit : public IFittingAlgorithm { public: /// Default constructor - Fit() : API::Algorithm(), m_domainType(API::IDomainCreator::Simple){}; + Fit() : IFittingAlgorithm(){} /// Algorithm's name for identification overriding a virtual method virtual const std::string name() const { return "Fit"; } /// Summary of algorithms purpose @@ -106,33 +107,11 @@ class DLLExport Fit : public API::Algorithm { /// Algorithm's version for identification overriding a virtual method virtual int version() const { return (1); } - /// Algorithm's category for identification overriding a virtual method - virtual const std::string category() const { return "Optimization"; } protected: - // Overridden Algorithm methods - void init(); - void exec(); - /// Override this method to perform a custom action right after a property was - /// set. - /// The argument is the property name. Default - do nothing. - virtual void afterPropertySet(const std::string &propName); - void setFunction(); - void addWorkspace(const std::string &workspaceNameProperty, - bool addProperties = true); - void addWorkspaces(); - /// Read domain type property and cache the value - void setDomainType(); + void initConcrete(); + void execConcrete(); void copyMinimizerOutput(const API::IFuncMinimizer &minimizer); - - /// Pointer to the fitting function - API::IFunction_sptr m_function; - /// Pointer to a domain creator - boost::shared_ptr m_domainCreator; - friend class API::IDomainCreator; - std::vector m_workspacePropertyNames; - /// Keep the domain type - API::IDomainCreator::DomainType m_domainType; }; } // namespace CurveFitting diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/IFittingAlgorithm.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/IFittingAlgorithm.h new file mode 100644 index 000000000000..cc3c026e3821 --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/IFittingAlgorithm.h @@ -0,0 +1,89 @@ +#ifndef MANTID_CURVEFITTING_IFITTINGALGORITHM_H_ +#define MANTID_CURVEFITTING_IFITTINGALGORITHM_H_ + +#include "MantidKernel/System.h" +#include "MantidAPI/Algorithm.h" +#include "MantidAPI/IDomainCreator.h" + +namespace Mantid { + +namespace API { + class IFunction; +} + +namespace CurveFitting { + +/** + + A base class for fitting algorithms. It declares two properties: + "Function" and "InputWorkspace". When these properties are set it + creates an appropriate domain creator. + + The concrete classes must implement methods: + - name() + - version() + - summary() + - exec() + - initConcrete() to declare more properties + + Copyright © 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge + National Laboratory & European Spallation Source + + 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: +*/ +class DLLExport IFittingAlgorithm : public API::Algorithm { +public: + IFittingAlgorithm(); + virtual ~IFittingAlgorithm(); + + virtual const std::string category() const; + +private: + void init(); + void exec(); + /// Child classes declare their properties here. + virtual void initConcrete() = 0; + /// Child classes implement the algorithm logic here. + virtual void execConcrete() = 0; + + virtual void afterPropertySet(const std::string &propName); + void addWorkspace(const std::string &workspaceNameProperty, + bool addProperties = true); + /// Read domain type property and cache the value + void setDomainType(); + +protected: + void setFunction(); + void addWorkspaces(); + + /// Keep the domain type + API::IDomainCreator::DomainType m_domainType; + /// Pointer to the fitting function + boost::shared_ptr m_function; + /// Pointer to a domain creator + boost::shared_ptr m_domainCreator; + std::vector m_workspacePropertyNames; + + friend class API::IDomainCreator; +}; + +} // namespace CurveFitting +} // namespace Mantid + +#endif /* MANTID_CURVEFITTING_IFITTINGALGORITHM_H_ */ \ No newline at end of file diff --git a/Code/Mantid/Framework/CurveFitting/src/CalculateChiSquared.cpp b/Code/Mantid/Framework/CurveFitting/src/CalculateChiSquared.cpp new file mode 100644 index 000000000000..f0fac91ba88d --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/src/CalculateChiSquared.cpp @@ -0,0 +1,115 @@ +#include "MantidCurveFitting/CalculateChiSquared.h" + +namespace Mantid { +namespace CurveFitting { + +using namespace Mantid::Kernel; +using namespace Mantid::API; + +// Register the algorithm into the AlgorithmFactory +DECLARE_ALGORITHM(CalculateChiSquared) + +//---------------------------------------------------------------------------------------------- +/** Constructor + */ +CalculateChiSquared::CalculateChiSquared() {} + +//---------------------------------------------------------------------------------------------- +/** Destructor + */ +CalculateChiSquared::~CalculateChiSquared() {} + +//---------------------------------------------------------------------------------------------- + +/// Algorithms name for identification. @see Algorithm::name +const std::string CalculateChiSquared::name() const { + return "CalculateChiSquared"; +} + +/// Algorithm's version for identification. @see Algorithm::version +int CalculateChiSquared::version() const { return 1; } + +/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary +const std::string CalculateChiSquared::summary() const { + return "Calculate chi squared for a function and a data set in a workspace."; +} + +//---------------------------------------------------------------------------------------------- +/// Initialize the algorithm's properties. +void CalculateChiSquared::initConcrete() { + declareProperty("ChiSquared", 0.0, "Output value of chi squared.", Direction::Output); + declareProperty("ChiSquaredDividedByDOF", 0.0, + "Output value of chi squared divided by the " + "number of degrees of freedom (NofData " + "- nOfParams).", Direction::Output); + declareProperty("ChiSquaredWeighted", 0.0, "Output value of weighted chi squared.", Direction::Output); + declareProperty("ChiSquaredWeightedDividedByDOF", 0.0, + "Output value of weighted chi squared divided by the " + "number of degrees of freedom (NofData " + "- nOfParams).", Direction::Output); +} + +//---------------------------------------------------------------------------------------------- +/// Execute the algorithm. +void CalculateChiSquared::execConcrete() { + + // Function may need some preparation. + m_function->setUpForFit(); + + API::FunctionDomain_sptr domain; + API::FunctionValues_sptr values; + m_domainCreator->createDomain(domain, values); + + // Do something with the function which may depend on workspace. + m_domainCreator->initFunction(m_function); + + // Calculate function values. + m_function->function(*domain, *values); + + // Get the number of free fitting parameters + size_t nParams = 0; + for (size_t i = 0; i < m_function->nParams(); ++i) { + if (!m_function->isFixed(i)) + nParams += 1; + } + double dof = - static_cast(nParams); + + // Calculate the chi squared. + double chiSquared = 0.0; + double chiSquaredWeighted = 0.0; + for(size_t i = 0; i < values->size(); ++i) { + auto weight = values->getFitWeight(i); + if (weight > 0.0) { + double tmp = values->getFitData(i) - values->getCalculated(i); + chiSquared += tmp * tmp; + tmp *= weight; + chiSquaredWeighted += tmp * tmp; + dof += 1.0; + } + } + g_log.notice() << "Chi squared " << chiSquared << std::endl; + g_log.notice() << "Chi squared weighted " << chiSquaredWeighted << std::endl; + + // Store the result. + setProperty("ChiSquared", chiSquared); + setProperty("chiSquaredWeighted", chiSquaredWeighted); + + // Divide by the DOF + g_log.debug() << "DOF " << dof << std::endl; + if (dof <= 0.0) { + dof = 1.0; + g_log.warning() << "DOF has a non-positive value, changing to 1.0." + << std::endl; + } + chiSquared /= dof; + chiSquaredWeighted /= dof; + g_log.notice() << "Chi squared / DOF " << chiSquared << std::endl; + g_log.notice() << "Chi squared weighed / DOF " << chiSquaredWeighted << std::endl; + + // Store the result. + setProperty("ChiSquaredDividedByDOF", chiSquared); + setProperty("ChiSquaredWeightedDividedByDOF", chiSquaredWeighted); +} + +} // namespace CurveFitting +} // namespace Mantid diff --git a/Code/Mantid/Framework/CurveFitting/src/Fit.cpp b/Code/Mantid/Framework/CurveFitting/src/Fit.cpp index 989d65bd9f8a..11ad2306a568 100644 --- a/Code/Mantid/Framework/CurveFitting/src/Fit.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/Fit.cpp @@ -2,312 +2,28 @@ // Includes //---------------------------------------------------------------------- #include "MantidCurveFitting/Fit.h" -#include "MantidCurveFitting/BoundaryConstraint.h" #include "MantidCurveFitting/CostFuncFitting.h" -#include "MantidCurveFitting/FitMW.h" -#include "MantidCurveFitting/MultiDomainCreator.h" -#include "MantidCurveFitting/Convolution.h" -#include "MantidCurveFitting/SeqDomainSpectrumCreator.h" -#include "MantidCurveFitting/LatticeDomainCreator.h" +#include "MantidAPI/CostFunctionFactory.h" #include "MantidAPI/FuncMinimizerFactory.h" +#include "MantidAPI/FunctionValues.h" #include "MantidAPI/IFuncMinimizer.h" -#include "MantidAPI/DomainCreatorFactory.h" -#include "MantidAPI/WorkspaceFactory.h" #include "MantidAPI/MatrixWorkspace.h" -#include "MantidAPI/IMDWorkspace.h" -#include "MantidAPI/FunctionProperty.h" -#include "MantidAPI/CostFunctionFactory.h" -#include "MantidAPI/CompositeDomain.h" -#include "MantidAPI/FunctionValues.h" -#include "MantidAPI/IFunctionMW.h" -#include "MantidAPI/IFunctionMD.h" -#include "MantidAPI/IFunction1DSpectrum.h" -#include "MantidAPI/ILatticeFunction.h" -#include "MantidAPI/MultiDomainFunction.h" -#include "MantidAPI/ITableWorkspace.h" - #include "MantidAPI/TableRow.h" -#include "MantidAPI/TextAxis.h" -#include "MantidKernel/UnitFactory.h" -#include "MantidKernel/Matrix.h" +#include "MantidAPI/WorkspaceFactory.h" + #include "MantidKernel/BoundedValidator.h" -#include "MantidKernel/ListValidator.h" #include "MantidKernel/StartsWithValidator.h" -#include -#include -#include - namespace Mantid { namespace CurveFitting { // Register the class into the algorithm factory DECLARE_ALGORITHM(Fit) -using API::IDomainCreator; - -/** - * Examine "Function" and "InputWorkspace" properties to decide which domain - * creator to use. - * @param propName :: A property name. - */ -void Fit::afterPropertySet(const std::string &propName) { - if (propName == "Function") { - setFunction(); - } else if (propName.size() >= 14 && - propName.substr(0, 14) == "InputWorkspace") { - if (getPointerToProperty("Function")->isDefault()) { - throw std::invalid_argument("Function must be set before InputWorkspace"); - } - addWorkspace(propName); - } else if (propName == "DomainType") { - setDomainType(); - } -} - -/** - * Read domain type property and cache the value - */ -void Fit::setDomainType() { - std::string domainType = getPropertyValue("DomainType"); - if (domainType == "Simple") { - m_domainType = IDomainCreator::Simple; - } else if (domainType == "Sequential") { - m_domainType = IDomainCreator::Sequential; - } else if (domainType == "Parallel") { - m_domainType = IDomainCreator::Parallel; - } else { - m_domainType = IDomainCreator::Simple; - } - Kernel::Property *prop = getPointerToProperty("Minimizer"); - auto minimizerProperty = - dynamic_cast *>(prop); - if (!minimizerProperty) { - throw std::logic_error("Failed to cast Property to PropertyWithValue"); - } - std::vector minimizerOptions = - API::FuncMinimizerFactory::Instance().getKeys(); - if (m_domainType != IDomainCreator::Simple) { - auto it = std::find(minimizerOptions.begin(), minimizerOptions.end(), - "Levenberg-Marquardt"); - minimizerOptions.erase(it); - } - minimizerProperty->replaceValidator(Kernel::IValidator_sptr( - new Kernel::StartsWithValidator(minimizerOptions))); -} - -/** - * Copy all output workspace properties from the minimizer to Fit algorithm. - * @param minimizer :: The minimizer to copy from. - */ -void Fit::copyMinimizerOutput(const API::IFuncMinimizer &minimizer) { - auto &properties = minimizer.getProperties(); - for (auto prop = properties.begin(); prop != properties.end(); ++prop) { - if ((**prop).direction() == Kernel::Direction::Output && - (**prop).isValid() == "") { - Kernel::Property *property = (**prop).clone(); - declareProperty(property); - } - } -} - -void Fit::setFunction() { - // get the function - m_function = getProperty("Function"); - auto mdf = boost::dynamic_pointer_cast(m_function); - if (mdf) { - size_t ndom = mdf->getMaxIndex() + 1; - m_workspacePropertyNames.resize(ndom); - m_workspacePropertyNames[0] = "InputWorkspace"; - for (size_t i = 1; i < ndom; ++i) { - std::string workspacePropertyName = - "InputWorkspace_" + boost::lexical_cast(i); - m_workspacePropertyNames[i] = workspacePropertyName; - if (!existsProperty(workspacePropertyName)) { - declareProperty( - new API::WorkspaceProperty( - workspacePropertyName, "", Kernel::Direction::Input), - "Name of the input Workspace"); - } - } - } else { - m_workspacePropertyNames.resize(1, "InputWorkspace"); - } -} - -/** - * Add a new workspace to the fit. The workspace is in the property named - * workspacePropertyName - * @param workspacePropertyName :: A workspace property name (eg InputWorkspace - * or InputWorkspace_2). - * The property must already exist in the algorithm. - * @param addProperties :: allow for declaration of properties that specify the - * dataset - * within the workspace to fit to. - */ -void Fit::addWorkspace(const std::string &workspacePropertyName, - bool addProperties) { - // get the workspace - API::Workspace_const_sptr ws = getProperty(workspacePropertyName); - // m_function->setWorkspace(ws); - const size_t n = std::string("InputWorkspace").size(); - const std::string suffix = - (workspacePropertyName.size() > n) ? workspacePropertyName.substr(n) : ""; - const size_t index = - suffix.empty() ? 0 : boost::lexical_cast(suffix.substr(1)); - - API::IFunction_sptr fun = getProperty("Function"); - IDomainCreator *creator = NULL; - setDomainType(); - - // ILatticeFunction requires API::LatticeDomain. - if (boost::dynamic_pointer_cast(fun)) { - creator = new LatticeDomainCreator(this, workspacePropertyName); - } else { - if (boost::dynamic_pointer_cast(ws) && - !boost::dynamic_pointer_cast(fun)) { - /* IFunction1DSpectrum needs a different domain creator. If a function - * implements that type, we need to react appropriately at this point. - * Otherwise, the default creator FitMW is used. - */ - if (boost::dynamic_pointer_cast(fun)) { - creator = new SeqDomainSpectrumCreator(this, workspacePropertyName); - } else { - creator = new FitMW(this, workspacePropertyName, m_domainType); - } - } else { - try { - creator = API::DomainCreatorFactory::Instance().createDomainCreator( - "FitMD", this, workspacePropertyName, m_domainType); - } - catch (Kernel::Exception::NotFoundError &) { - throw std::invalid_argument("Unsupported workspace type" + ws->id()); - } - } - } - - if (!m_domainCreator) { - if (m_workspacePropertyNames.empty()) { - // this defines the function and fills in m_workspacePropertyNames with - // names of the sort InputWorkspace_# - setFunction(); - } - auto multiFun = boost::dynamic_pointer_cast(fun); - if (multiFun) { - auto multiCreator = - new MultiDomainCreator(this, m_workspacePropertyNames); - multiCreator->setCreator(index, creator); - m_domainCreator.reset(multiCreator); - creator->declareDatasetProperties(suffix, addProperties); - } else { - m_domainCreator.reset(creator); - creator->declareDatasetProperties(suffix, addProperties); - } - } else { - boost::shared_ptr multiCreator = - boost::dynamic_pointer_cast(m_domainCreator); - if (!multiCreator) { - throw std::runtime_error( - std::string("MultiDomainCreator expected, found ") + - typeid(*m_domainCreator.get()).name()); - } - if (!multiCreator->hasCreator(index)) { - creator->declareDatasetProperties(suffix, addProperties); - } - multiCreator->setCreator(index, creator); - } -} - -/** - * Collect all input workspace property names in the m_workspacePropertyNames - * vector - */ -void Fit::addWorkspaces() { - setDomainType(); - auto multiFun = - boost::dynamic_pointer_cast(m_function); - if (multiFun) { - m_domainCreator.reset( - new MultiDomainCreator(this, m_workspacePropertyNames)); - } - auto props = getProperties(); - for (auto prop = props.begin(); prop != props.end(); ++prop) { - if ((**prop).direction() == Kernel::Direction::Input && - dynamic_cast(*prop)) { - const std::string workspacePropertyName = (**prop).name(); - API::Workspace_const_sptr ws = getProperty(workspacePropertyName); - IDomainCreator *creator = NULL; - - // ILatticeFunction requires API::LatticeDomain. - if (boost::dynamic_pointer_cast(m_function)) { - creator = new LatticeDomainCreator(this, workspacePropertyName); - } else { - if (boost::dynamic_pointer_cast(ws) && - !boost::dynamic_pointer_cast(m_function)) { - /* IFunction1DSpectrum needs a different domain creator. If a function - * implements that type, we need to react appropriately at this point. - * Otherwise, the default creator FitMW is used. - */ - if (boost::dynamic_pointer_cast( - m_function)) { - creator = new SeqDomainSpectrumCreator(this, workspacePropertyName); - } else { - creator = new FitMW(this, workspacePropertyName, m_domainType); - } - } else { // don't know what to do with this workspace - try { - creator = API::DomainCreatorFactory::Instance().createDomainCreator( - "FitMD", this, workspacePropertyName, m_domainType); - } - catch (Kernel::Exception::NotFoundError &) { - throw std::invalid_argument("Unsupported workspace type" + - ws->id()); - } - } - } - - const size_t n = std::string("InputWorkspace").size(); - const std::string suffix = (workspacePropertyName.size() > n) - ? workspacePropertyName.substr(n) - : ""; - const size_t index = - suffix.empty() ? 0 : boost::lexical_cast(suffix.substr(1)); - creator->declareDatasetProperties(suffix, false); - m_workspacePropertyNames.push_back(workspacePropertyName); - if (!m_domainCreator) { - m_domainCreator.reset(creator); - } - auto multiCreator = - boost::dynamic_pointer_cast(m_domainCreator); - if (multiCreator) { - multiCreator->setCreator(index, creator); - } - } - } -} - /** Initialisation method */ -void Fit::init() { - declareProperty( - new API::FunctionProperty("Function"), - "Parameters defining the fitting function and its initial values"); - - declareProperty(new API::WorkspaceProperty( - "InputWorkspace", "", Kernel::Direction::Input), - "Name of the input Workspace"); - - std::vector domainTypes; - domainTypes.push_back("Simple"); - domainTypes.push_back("Sequential"); - domainTypes.push_back("Parallel"); - declareProperty( - "DomainType", "Simple", - Kernel::IValidator_sptr( - new Kernel::ListValidator(domainTypes)), - "The type of function domain to use: Simple, Sequential, or Parallel.", - Kernel::Direction::Input); +void Fit::initConcrete() { declareProperty("Ties", "", Kernel::Direction::Input); getPointerToProperty("Ties") @@ -321,8 +37,6 @@ void Fit::init() { declareProperty( "MaxIterations", 500, mustBePositive->clone(), "Stop after this number of iterations if a good fit is not found"); - declareProperty("IgnoreInvalidData", false, - "Flag to ignore infinities, NaNs and data with zero errors."); declareProperty("OutputStatus", "", Kernel::Direction::Output); getPointerToProperty("OutputStatus") ->setDocumentation("Whether the fit was successful"); @@ -383,16 +97,26 @@ void Fit::init() { "Output is an empty string)."); } +/** + * Copy all output workspace properties from the minimizer to Fit algorithm. + * @param minimizer :: The minimizer to copy from. + */ +void Fit::copyMinimizerOutput(const API::IFuncMinimizer &minimizer) { + auto &properties = minimizer.getProperties(); + for (auto prop = properties.begin(); prop != properties.end(); ++prop) { + if ((**prop).direction() == Kernel::Direction::Output && + (**prop).isValid() == "") { + Kernel::Property *property = (**prop).clone(); + declareProperty(property); + } + } +} + /** Executes the algorithm * * @throw runtime_error Thrown if algorithm cannot execute */ -void Fit::exec() { - // this is to make it work with AlgorithmProxy - if (!m_domainCreator) { - setFunction(); - addWorkspaces(); - } +void Fit::execConcrete() { std::string ties = getPropertyValue("Ties"); if (!ties.empty()) { @@ -407,8 +131,7 @@ void Fit::exec() { m_function->setUpForFit(); API::FunctionDomain_sptr domain; - API::FunctionValues_sptr values; // TODO: should values be part of domain? - m_domainCreator->ignoreInvalidData(getProperty("IgnoreInvalidData")); + API::FunctionValues_sptr values; m_domainCreator->createDomain(domain, values); // do something with the function which may depend on workspace diff --git a/Code/Mantid/Framework/CurveFitting/src/IFittingAlgorithm.cpp b/Code/Mantid/Framework/CurveFitting/src/IFittingAlgorithm.cpp new file mode 100644 index 000000000000..efb055c2aabc --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/src/IFittingAlgorithm.cpp @@ -0,0 +1,289 @@ +#include "MantidCurveFitting/IFittingAlgorithm.h" + +#include "MantidCurveFitting/FitMW.h" +#include "MantidCurveFitting/MultiDomainCreator.h" +#include "MantidCurveFitting/SeqDomainSpectrumCreator.h" +#include "MantidCurveFitting/LatticeDomainCreator.h" + +#include "MantidAPI/FunctionProperty.h" +#include "MantidAPI/IFunctionMD.h" +#include "MantidAPI/IFunction1DSpectrum.h" +#include "MantidAPI/ILatticeFunction.h" +#include "MantidAPI/MultiDomainFunction.h" + +#include "MantidKernel/ListValidator.h" + + +namespace Mantid { +namespace CurveFitting { + +using namespace Mantid::Kernel; +using namespace Mantid::API; + +//---------------------------------------------------------------------------------------------- +namespace { +/// Create a domain creator for a particular function and workspace pair. +IDomainCreator *createDomainCreator(const IFunction *fun, const Workspace *ws, + const std::string &workspacePropertyName, + IPropertyManager* manager, + IDomainCreator::DomainType domainType + ) { + + IDomainCreator *creator = NULL; + + // ILatticeFunction requires API::LatticeDomain. + if (dynamic_cast(fun)) { + creator = new LatticeDomainCreator(manager, workspacePropertyName); + } else { + if (dynamic_cast(ws) && + !dynamic_cast(fun)) { + /* IFunction1DSpectrum needs a different domain creator. If a function + * implements that type, we need to react appropriately at this point. + * Otherwise, the default creator FitMW is used. + */ + if (dynamic_cast(fun)) { + creator = new SeqDomainSpectrumCreator(manager, workspacePropertyName); + } else { + creator = new FitMW(manager, workspacePropertyName, domainType); + } + } else { + try { + creator = API::DomainCreatorFactory::Instance().createDomainCreator( + "FitMD", manager, workspacePropertyName, domainType); + } + catch (Kernel::Exception::NotFoundError &) { + throw std::invalid_argument("Unsupported workspace type" + ws->id()); + } + } + } + return creator; +} +} + +//---------------------------------------------------------------------------------------------- +/** Constructor + */ +IFittingAlgorithm::IFittingAlgorithm() + : API::Algorithm(), m_domainType(API::IDomainCreator::Simple) {} + +//---------------------------------------------------------------------------------------------- +/** Destructor + */ +IFittingAlgorithm::~IFittingAlgorithm() {} + +/// Algorithm's category for identification. @see Algorithm::category +const std::string IFittingAlgorithm::category() const { return "Optimization"; } + +//---------------------------------------------------------------------------------------------- +/** Initialize the algorithm's properties. + */ +void IFittingAlgorithm::init() { + declareProperty( + new API::FunctionProperty("Function"), + "Parameters defining the fitting function and its initial values"); + + declareProperty(new API::WorkspaceProperty( + "InputWorkspace", "", Kernel::Direction::Input), + "Name of the input Workspace"); + declareProperty("IgnoreInvalidData", false, + "Flag to ignore infinities, NaNs and data with zero errors."); + + std::vector domainTypes; + domainTypes.push_back("Simple"); + domainTypes.push_back("Sequential"); + domainTypes.push_back("Parallel"); + declareProperty( + "DomainType", "Simple", + Kernel::IValidator_sptr( + new Kernel::ListValidator(domainTypes)), + "The type of function domain to use: Simple, Sequential, or Parallel.", + Kernel::Direction::Input); + + initConcrete(); +} + +/** + * Examine "Function" and "InputWorkspace" properties to decide which domain + * creator to use. + * @param propName :: A property name. + */ +void IFittingAlgorithm::afterPropertySet(const std::string &propName) { + if (propName == "Function") { + setFunction(); + } else if (propName.size() >= 14 && + propName.substr(0, 14) == "InputWorkspace") { + if (getPointerToProperty("Function")->isDefault()) { + throw std::invalid_argument("Function must be set before InputWorkspace"); + } + addWorkspace(propName); + } else if (propName == "DomainType") { + setDomainType(); + } +} + +/** + * Read domain type property and cache the value + */ +void IFittingAlgorithm::setDomainType() { + std::string domainType = getPropertyValue("DomainType"); + if (domainType == "Simple") { + m_domainType = IDomainCreator::Simple; + } else if (domainType == "Sequential") { + m_domainType = IDomainCreator::Sequential; + } else if (domainType == "Parallel") { + m_domainType = IDomainCreator::Parallel; + } else { + m_domainType = IDomainCreator::Simple; + } + //Kernel::Property *prop = getPointerToProperty("Minimizer"); + //auto minimizerProperty = + // dynamic_cast *>(prop); + //std::vector minimizerOptions = + // API::FuncMinimizerFactory::Instance().getKeys(); + //if (m_domainType != IDomainCreator::Simple) { + // auto it = std::find(minimizerOptions.begin(), minimizerOptions.end(), + // "Levenberg-Marquardt"); + // minimizerOptions.erase(it); + //} + //minimizerProperty->replaceValidator(Kernel::IValidator_sptr( + // new Kernel::StartsWithValidator(minimizerOptions))); +} + +void IFittingAlgorithm::setFunction() { + // get the function + m_function = getProperty("Function"); + auto mdf = boost::dynamic_pointer_cast(m_function); + if (mdf) { + size_t ndom = mdf->getMaxIndex() + 1; + m_workspacePropertyNames.resize(ndom); + m_workspacePropertyNames[0] = "InputWorkspace"; + for (size_t i = 1; i < ndom; ++i) { + std::string workspacePropertyName = + "InputWorkspace_" + boost::lexical_cast(i); + m_workspacePropertyNames[i] = workspacePropertyName; + if (!existsProperty(workspacePropertyName)) { + declareProperty( + new API::WorkspaceProperty( + workspacePropertyName, "", Kernel::Direction::Input), + "Name of the input Workspace"); + } + } + } else { + m_workspacePropertyNames.resize(1, "InputWorkspace"); + } +} + +/** + * Add a new workspace to the fit. The workspace is in the property named + * workspacePropertyName. + * @param workspacePropertyName :: A workspace property name (eg InputWorkspace + * or InputWorkspace_2). The property must already exist in the algorithm. + * @param addProperties :: allow for declaration of properties that specify the + * dataset within the workspace to fit to. + */ +void IFittingAlgorithm::addWorkspace(const std::string &workspacePropertyName, + bool addProperties) { + // get the workspace + API::Workspace_const_sptr ws = getProperty(workspacePropertyName); + // m_function->setWorkspace(ws); + const size_t n = std::string("InputWorkspace").size(); + const std::string suffix = + (workspacePropertyName.size() > n) ? workspacePropertyName.substr(n) : ""; + const size_t index = + suffix.empty() ? 0 : boost::lexical_cast(suffix.substr(1)); + + IFunction_sptr fun = getProperty("Function"); + setDomainType(); + + IDomainCreator *creator = createDomainCreator( + fun.get(), ws.get(), workspacePropertyName, this, m_domainType); + + if (!m_domainCreator) { + if (m_workspacePropertyNames.empty()) { + // this defines the function and fills in m_workspacePropertyNames with + // names of the sort InputWorkspace_# + setFunction(); + } + auto multiFun = boost::dynamic_pointer_cast(fun); + if (multiFun) { + auto multiCreator = + new MultiDomainCreator(this, m_workspacePropertyNames); + multiCreator->setCreator(index, creator); + m_domainCreator.reset(multiCreator); + creator->declareDatasetProperties(suffix, addProperties); + } else { + m_domainCreator.reset(creator); + creator->declareDatasetProperties(suffix, addProperties); + } + } else { + boost::shared_ptr multiCreator = + boost::dynamic_pointer_cast(m_domainCreator); + if (!multiCreator) { + throw std::runtime_error( + std::string("MultiDomainCreator expected, found ") + + typeid(*m_domainCreator.get()).name()); + } + if (!multiCreator->hasCreator(index)) { + creator->declareDatasetProperties(suffix, addProperties); + } + multiCreator->setCreator(index, creator); + } +} + +/** + * Collect all input workspace property names in the m_workspacePropertyNames + * vector. + */ +void IFittingAlgorithm::addWorkspaces() { + setDomainType(); + auto multiFun = + boost::dynamic_pointer_cast(m_function); + if (multiFun) { + m_domainCreator.reset( + new MultiDomainCreator(this, m_workspacePropertyNames)); + } + auto props = getProperties(); + for (auto prop = props.begin(); prop != props.end(); ++prop) { + if ((**prop).direction() == Kernel::Direction::Input && + dynamic_cast(*prop)) { + const std::string workspacePropertyName = (**prop).name(); + API::Workspace_const_sptr ws = getProperty(workspacePropertyName); + IDomainCreator *creator = + createDomainCreator(m_function.get(), ws.get(), workspacePropertyName, + this, m_domainType); + + const size_t n = std::string("InputWorkspace").size(); + const std::string suffix = (workspacePropertyName.size() > n) + ? workspacePropertyName.substr(n) + : ""; + const size_t index = + suffix.empty() ? 0 : boost::lexical_cast(suffix.substr(1)); + creator->declareDatasetProperties(suffix, false); + m_workspacePropertyNames.push_back(workspacePropertyName); + if (!m_domainCreator) { + m_domainCreator.reset(creator); + } + auto multiCreator = + boost::dynamic_pointer_cast(m_domainCreator); + if (multiCreator) { + multiCreator->setCreator(index, creator); + } + } + } +} +//---------------------------------------------------------------------------------------------- +/// Execute the algorithm. +void IFittingAlgorithm::exec() { + + // This is to make it work with AlgorithmProxy + if (!m_domainCreator) { + setFunction(); + addWorkspaces(); + } + m_domainCreator->ignoreInvalidData(getProperty("IgnoreInvalidData")); + // Execute the concrete algorithm. + this->execConcrete(); +} + +} // namespace CurveFitting +} // namespace Mantid \ No newline at end of file diff --git a/Code/Mantid/Framework/CurveFitting/test/CalculateChiSquaredTest.h b/Code/Mantid/Framework/CurveFitting/test/CalculateChiSquaredTest.h new file mode 100644 index 000000000000..0f428ac39a87 --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/test/CalculateChiSquaredTest.h @@ -0,0 +1,362 @@ +#ifndef MANTID_CURVEFITTING_CALCULATECHISQUAREDTEST_H_ +#define MANTID_CURVEFITTING_CALCULATECHISQUAREDTEST_H_ + +#include + +#include "MantidCurveFitting/CalculateChiSquared.h" +#include "MantidAPI/FunctionDomain1D.h" +#include "MantidAPI/FunctionFactory.h" +#include "MantidAPI/FunctionValues.h" +#include "MantidAPI/IFunction1D.h" +#include "MantidKernel/EmptyValues.h" + +#include +#include +#include + +using Mantid::CurveFitting::CalculateChiSquared; +using namespace Mantid; +using namespace Mantid::API; + +class CalculateChiSquaredTest : 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 CalculateChiSquaredTest *createSuite() { + return new CalculateChiSquaredTest(); + } + static void destroySuite(CalculateChiSquaredTest *suite) { delete suite; } + + void test_Init() { + CalculateChiSquared alg; + TS_ASSERT_THROWS_NOTHING(alg.initialize()) + TS_ASSERT(alg.isInitialized()) + } + + void test_1D_empty_defaults() { + Tester tester; + tester.set1DFunction(); + tester.set1DSpectrumEmpty(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 20338.0, 1.0); + } + + void test_1D_empty_all_x_range() { + Tester tester; + tester.set1DFunction(); + tester.set1DSpectrumEmpty(); + tester.setXRange_All(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 20338.0, 1.0); + } + + void test_1D_empty_smaller() { + Tester tester; + tester.set1DFunction(); + tester.set1DSpectrumEmpty(); + tester.setXRange_Smaller(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 1189.0, 1.0); + } + + void test_1D_empty_smaller1() { + Tester tester; + tester.set1DFunction(); + tester.set1DSpectrumEmpty(); + tester.setXRange_Smaller_bin_boundaries(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 1189.0, 1.0); + } + + void test_1D_values() { + Tester tester; + tester.set1DFunction(); + tester.set1DSpectrumValues(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 1655.0, 1.0); + } + + void test_1D_values_smaller() { + Tester tester; + tester.set1DFunction(); + tester.set1DSpectrumValues(); + tester.setXRange_Smaller(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 153.0, 1.0); + } + + void test_1D_values_point_data() { + Tester tester(3,10,false); + tester.set1DFunction(); + tester.set1DSpectrumValues(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 307.0, 1.0); + } + + void test_1D_workspace_index() { + Tester tester; + tester.set1DFunction(); + tester.set1DSpectrumValues(5); + tester.setXRange_Smaller(); + tester.setWorkspaceIndex(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 151.0, 1.0); + } + + void test_1D_values_dont_ignore_invalid() { + Tester tester; + tester.set1DFunction(); + tester.set1DSpectrumValuesInvalid(); + tester.runAlgorithm(); + tester.checkFailed(); + } + + void test_1D_values_ignore_invalid() { + Tester tester; + tester.set1DFunction(); + tester.set1DSpectrumValuesInvalid(); + tester.setIgnoreInvalidData(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 1450.39, 1.0); + } + + void test_1D_values_divide_by_dof_fixed_params() { + Tester tester(1); + tester.set1DFunction(); + tester.set1DSpectrumValues(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 1655.0, 1.0); + } + + void test_1D_values_divide_by_dof_zero() { + Tester tester(3,3); + tester.set1DFunction(); + tester.set1DSpectrumValues(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 5069.0, 1.0); + } + + void test_1D_values_divide_by_dof_negative() { + Tester tester(3,2); + tester.set1DFunction(); + tester.set1DSpectrumValues(); + tester.runAlgorithm(); + tester.check1DSpectrum(); + TS_ASSERT_DELTA(tester.chiSquared, 7151.0, 1.0); + } + +private: + class Tester { + // input parameters + size_t nParams; + size_t nData; + bool isHisto; + double xMin; + double xMax; + std::vector xBins; + std::vector xValues; + + // values for algorithm input properties + IFunction_sptr function; + Workspace_sptr workspace; + int workspaceIndex; + double StartX; + double EndX; + bool ignoreInvalidData; + + void makeXValues() { + size_t dlt = isHisto ? 1 : 0; + xBins.resize(nData + dlt); + double dx = (xMax - xMin) / double(xBins.size() - 1); + for (size_t i = 0; i < xBins.size(); ++i) { + xBins[i] = xMin + double(i) * dx; + } + if (isHisto) { + xValues.resize(nData); + std::transform(xBins.begin(), xBins.end() - 1, xValues.begin(), + std::bind2nd(std::plus(), dx / 2)); + } else { + xValues = xBins; + } + } + + void setDefaultXRange() { + if (StartX == EMPTY_DBL()) + StartX = xMin; + if (EndX == EMPTY_DBL()) { + EndX = xMax; + } else { + auto ix = std::upper_bound(xBins.begin(), xBins.end(), EndX); + if (ix != xBins.end()) EndX = *ix; + else + EndX = xMax; + } + } + + bool isGoodValue(double y, double e) { + return !ignoreInvalidData || + (!boost::math::isnan(y) && !boost::math::isinf(y) && + !boost::math::isnan(e) && !boost::math::isinf(e) && e > 0); + } + + public: + // algorithm output + double chiSquared; + double chiSquaredDividedByDOF; + double chiSquaredWeighted; + double chiSquaredWeightedDividedByDOF; + bool isExecuted; + + Tester(size_t np = 3, size_t nd = 10, bool histo = true) + : nParams(np), nData(nd), isHisto(histo), xMin(-10), xMax(10), + workspaceIndex(0), StartX(EMPTY_DBL()), EndX(EMPTY_DBL()), + ignoreInvalidData(false), chiSquared(-1), chiSquaredDividedByDOF(-1), + chiSquaredWeighted(-1), chiSquaredWeightedDividedByDOF(-1), + isExecuted(false) { + makeXValues(); + } + + void runAlgorithm() { + CalculateChiSquared alg; + TS_ASSERT_THROWS_NOTHING(alg.initialize()) + TS_ASSERT(alg.isInitialized()) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("Function", function)); + TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", workspace)); + TS_ASSERT_THROWS_NOTHING(alg.setProperty("IgnoreInvalidData", ignoreInvalidData)); + if (dynamic_cast(function.get())) { + TS_ASSERT_THROWS_NOTHING(alg.setProperty("WorkspaceIndex", workspaceIndex)); + TS_ASSERT_THROWS_NOTHING(alg.setProperty("StartX", StartX)); + TS_ASSERT_THROWS_NOTHING(alg.setProperty("EndX", EndX)); + } + TS_ASSERT_THROWS_NOTHING(alg.execute()); + isExecuted = alg.isExecuted(); + if (isExecuted) { + chiSquared = alg.getProperty("ChiSquared"); + chiSquaredDividedByDOF = alg.getProperty("ChiSquaredDividedByDOF"); + chiSquaredWeighted = alg.getProperty("ChiSquaredWeighted"); + chiSquaredWeightedDividedByDOF = alg.getProperty("ChiSquaredWeightedDividedByDOF"); + } + } + + void setXRange_All() { + StartX = xMin; + EndX = xMax; + } + + void setXRange_Smaller_bin_boundaries() { + StartX = xBins[3]; + EndX = xBins[7]; + } + + void setXRange_Smaller() { + StartX = xBins[3] - 0.3; + EndX = xBins[7] + 0.7; + } + + void setWorkspaceIndex() { + workspaceIndex = 3; + } + + void setIgnoreInvalidData() { + ignoreInvalidData = true; + } + + void set1DFunction() { + const std::string fun = + "name=UserFunction,Formula=a+b*x+c*x^2,a=1,b=1,c=1"; + function = FunctionFactory::Instance().createInitialized(fun); + if (nParams < function->nParams()) { + const size_t dn = function->nParams() - nParams; + for(size_t i = 0; i < dn; ++i) { + function->fix(i); + } + } else { + TS_ASSERT_EQUALS(function->nParams(), nParams); + } + } + + void set1DSpectrumEmpty() { + + const size_t nSpec = 1; + size_t dn = isHisto ? 1 : 0; + auto space = WorkspaceFactory::Instance().create("Workspace2D", nSpec, + nData + dn, nData); + space->dataX(0).assign(xBins.begin(), xBins.end()); + space->dataE(0).assign(nData, 10.0); + workspace = space; + } + + void set1DSpectrumValues(const size_t nSpec = 1) { + size_t dn = isHisto ? 1 : 0; + auto space = WorkspaceFactory::Instance().create("Workspace2D", nSpec, + nData + dn, nData); + for (size_t spec = 0; spec < nSpec; ++spec) { + space->dataX(spec).assign(xBins.begin(), xBins.end()); + for (size_t i = 0; i < nData; ++i) { + const double x = space->readX(0)[i]; + space->dataY(spec)[i] = (1.1 + 0.1 * double(spec)) * (1.0 + x + x * x); + space->dataE(spec)[i] = 10.0; + } + } + workspace = space; + } + + void set1DSpectrumValuesInvalid() { + set1DSpectrumValues(); + auto &yValues = dynamic_cast(*workspace).dataY(workspaceIndex); + yValues[2] = std::numeric_limits::infinity(); + yValues[4] = std::numeric_limits::quiet_NaN(); + auto &eValues = dynamic_cast(*workspace).dataE(workspaceIndex); + eValues[6] = -1; + } + + void check1DSpectrum() { + TS_ASSERT(isExecuted); + setDefaultXRange(); + double sum2 = 0.0; + double sum2w = 0.0; + auto &yValues = dynamic_cast(*workspace).readY(workspaceIndex); + auto &eValues = dynamic_cast(*workspace).readE(workspaceIndex); + double dof = - double(nParams); + for (size_t i = 0; i < xValues.size(); ++i) { + const double xValue = xValues[i]; + if (xValue >= StartX && xValue <= EndX && isGoodValue(yValues[i],eValues[i])) { + FunctionDomain1DVector x(xValue); + FunctionValues y(x); + function->function(x, y); + double tmp = yValues[i] - y[0]; + //std::cerr << "test " << xValue << ' ' << yValues[i] << ' ' << y[0] << std::endl; + sum2 += tmp * tmp; + tmp /= eValues[i]; + sum2w += tmp * tmp; + dof += 1.0; + } + } + TS_ASSERT_DIFFERS(sum2, 0); + TS_ASSERT_DELTA(sum2, chiSquared, 1e-10); + TS_ASSERT_DELTA(sum2w, chiSquaredWeighted, 1e-10); + if (dof <= 0.0) dof = 1.0; + sum2 /= dof; + sum2w /= dof; + TS_ASSERT_DELTA(sum2, chiSquaredDividedByDOF, 1e-10); + TS_ASSERT_DELTA(sum2w, chiSquaredWeightedDividedByDOF, 1e-10); + } + + void checkFailed() { + TS_ASSERT(!isExecuted); + } + }; +}; + +#endif /* MANTID_CURVEFITTING_CALCULATECHISQUAREDTEST_H_ */ diff --git a/Code/Mantid/Framework/PythonInterface/mantid/simpleapi.py b/Code/Mantid/Framework/PythonInterface/mantid/simpleapi.py index d2186a1eac62..ed32cd3acc74 100644 --- a/Code/Mantid/Framework/PythonInterface/mantid/simpleapi.py +++ b/Code/Mantid/Framework/PythonInterface/mantid/simpleapi.py @@ -32,7 +32,7 @@ #------------------------ Specialized function calls -------------------------- # List of specialized algorithms -__SPECIALIZED_FUNCTIONS__ = ["Load", "Fit", "CutMD"] +__SPECIALIZED_FUNCTIONS__ = ["Load", "CutMD"] # List of specialized algorithms __MDCOORD_FUNCTIONS__ = ["PeakIntensityVsRadius", "CentroidPeaksMD","IntegratePeaksMD"] # The "magic" keyword to enable/disable logging @@ -171,9 +171,71 @@ def LoadDialog(*args, **kwargs): #---------------------------- Fit --------------------------------------------- +def fitting_algorithm(f): + """ + Decorator generating code for fitting algorithms (currently Fit and CalculateChiSquared). + When applied to a function definition this decorator replaces its code with code of + function 'wrapper' defined below. + """ + + function_name = f.__name__ + + def wrapper(*args, **kwargs): + Function, InputWorkspace = _get_mandatory_args(function_name, ["Function", "InputWorkspace"], *args, **kwargs) + # Remove from keywords so it is not set twice + if "Function" in kwargs: + del kwargs['Function'] + if "InputWorkspace" in kwargs: + del kwargs['InputWorkspace'] + + # Check for behaviour consistent with old API + if type(Function) == str and Function in _api.AnalysisDataService: + raise ValueError("Fit API has changed. The function must now come first in the argument list and the workspace second.") + # Create and execute + algm = _create_algorithm_object(function_name) + _set_logging_option(algm, kwargs) + algm.setProperty('Function', Function) # Must be set first + algm.setProperty('InputWorkspace', InputWorkspace) + + # Set all workspace properties before others + for key in kwargs.keys(): + if key.startswith('InputWorkspace_'): + algm.setProperty(key, kwargs[key]) + del kwargs[key] + + lhs = _kernel.funcreturns.lhs_info() + # Check for any properties that aren't known and warn they will not be used + for key in kwargs.keys(): + if key not in algm: + logger.warning("You've passed a property (%s) to %s() that doesn't apply to any of the input workspaces." % (key,function_name)) + del kwargs[key] + set_properties(algm, **kwargs) + algm.execute() + + return _gather_returns(function_name, lhs, algm) + + # Have a better load signature for autocomplete + _signature = "\bFunction, InputWorkspace" + # Getting the code object for Load + _f = wrapper.func_code + # Creating a new code object nearly identical, but with the two variable names replaced + # by the property list. + _c = _f.__new__(_f.__class__, _f.co_argcount, _f.co_nlocals, _f.co_stacksize, _f.co_flags, _f.co_code, _f.co_consts, _f.co_names,\ + (_signature, "kwargs"), _f.co_filename, _f.co_name, _f.co_firstlineno, _f.co_lnotab, _f.co_freevars) + + # Replace the code object of the wrapper function + wrapper.func_code = _c + wrapper.__doc__ = f.__doc__ + if not function_name in __SPECIALIZED_FUNCTIONS__: + __SPECIALIZED_FUNCTIONS__.append(function_name) + + return wrapper + +# Use a python decorator (defined above) to generate the code for this function. +@fitting_algorithm def Fit(*args, **kwargs): """ - Fit defines the interface to the fitting 562 within Mantid. + Fit defines the interface to the fitting within Mantid. It can work with arbitrary data sources and therefore some options are only available when the function & workspace type are known. @@ -185,50 +247,21 @@ def Fit(*args, **kwargs): Fit(Function='name=LinearBackground,A0=0.3', InputWorkspace=dataWS', StartX='0.05',EndX='1.0',Output="Z1") """ - Function, InputWorkspace = _get_mandatory_args('Fit', ["Function", "InputWorkspace"], *args, **kwargs) - # Remove from keywords so it is not set twice - if "Function" in kwargs: - del kwargs['Function'] - if "InputWorkspace" in kwargs: - del kwargs['InputWorkspace'] + return None - # Check for behaviour consistent with old API - if type(Function) == str and Function in _api.AnalysisDataService: - raise ValueError("Fit API has changed. The function must now come first in the argument list and the workspace second.") - # Create and execute - algm = _create_algorithm_object('Fit') - _set_logging_option(algm, kwargs) - algm.setProperty('Function', Function) # Must be set first - algm.setProperty('InputWorkspace', InputWorkspace) - - # Set all workspace properties before others - for key in kwargs.keys(): - if key.startswith('InputWorkspace_'): - algm.setProperty(key, kwargs[key]) - del kwargs[key] - - lhs = _kernel.funcreturns.lhs_info() - # Check for any properties that aren't known and warn they will not be used - for key in kwargs.keys(): - if key not in algm: - logger.warning("You've passed a property (%s) to Fit() that doesn't apply to any of the input workspaces." % key) - del kwargs[key] - set_properties(algm, **kwargs) - algm.execute() - - return _gather_returns('Fit', lhs, algm) - -# Have a better load signature for autocomplete -_signature = "\bFunction, InputWorkspace" -# Getting the code object for Load -_f = Fit.func_code -# Creating a new code object nearly identical, but with the two variable names replaced -# by the property list. -_c = _f.__new__(_f.__class__, _f.co_argcount, _f.co_nlocals, _f.co_stacksize, _f.co_flags, _f.co_code, _f.co_consts, _f.co_names,\ - (_signature, "kwargs"), _f.co_filename, _f.co_name, _f.co_firstlineno, _f.co_lnotab, _f.co_freevars) +# Use a python decorator (defined above) to generate the code for this function. +@fitting_algorithm +def CalculateChiSquared(*args, **kwargs): + """ + This function calculates chi squared claculation for a function and a data set. + The data set is defined in a way similar to Fit algorithm. -# Replace the code object of the wrapper function -Fit.func_code = _c + Example: + chi2_1, chi2_2, chi2_3, chi2_4 = \\ + CalculateChiSquared(Function='name=LinearBackground,A0=0.3', InputWorkspace=dataWS', + StartX='0.05',EndX='1.0') + """ + return None def FitDialog(*args, **kwargs): """Popup a dialog for the Load algorithm. More help on the Load function diff --git a/Code/Mantid/docs/source/algorithms/CalculateChiSquared-v1.rst b/Code/Mantid/docs/source/algorithms/CalculateChiSquared-v1.rst new file mode 100644 index 000000000000..6e361be09935 --- /dev/null +++ b/Code/Mantid/docs/source/algorithms/CalculateChiSquared-v1.rst @@ -0,0 +1,88 @@ + +.. algorithm:: + +.. summary:: + +.. alias:: + +.. properties:: + +Description +----------- + +Calculates a measure of the goodness of a fit in four ways. + +The ChiSquared property returns the sum of squares of differences between the calculated and measured values: + +:math:`\chi_{1}^{2} = \sum_{i} (y_i - f_i)^2` + +where :math:`y_i` and :math:`f_i` are the measured and calculated values at i-th point. + +The ChiSquaredDividedByDOF is ChiSquared divided by the number of degrees of freedom (DOF): + +:math:`\chi_{2}^{2} = \frac{1}{DOF}\sum_{i} (y_i - f_i)^2` + +:math:`DOF = N_d - N_p` where :math:`N_d` is the number of data points used in fitting and :math:`N_p` +is the number of free (not fixed or tied) parameters of the function. + +The ChiSquaredWeighted property sums the squares of the differences divided by the data errors: + +:math:`\chi_{3}^{2} = \sum_{i} \left(\frac{y_i - f_i}{\sigma_i}\right)^2` + +Finally, ChiSquaredWeightedDividedByDOF is + +:math:`\chi_{4}^{2} = \chi_{3}^{2} / DOF` + +Usage +----- +**Example - CalculateChiSquared** + +.. testcode:: CalculateChiSquaredExample + + import numpy as np + + # Create a data set + x = np.linspace(0,1,10) + y = 1.0 + 2.0 * x + e = np.sqrt(y) + ws = CreateWorkspace(DataX=x, DataY=y, DataE=e) + + # Define a function + func = 'name=LinearBackground,A0=1.1,A1=1.9' + + # Calculate the chi squared + chi2,chi2dof,chi2W,chi2Wdof = CalculateChiSquared(func,ws) + + print 'Chi squared is %s' % chi2 + print 'Chi squared / DOF is %s' % chi2dof + print 'Chi squared weighted is %s' % chi2W + print 'Chi squared weighted / DOF is %s' % chi2Wdof + print + + # Define a function that models the data exactly + func = 'name=LinearBackground,A0=1.0,A1=2.0' + + # Calculate the chi squared + chi2,chi2dof,chi2W,chi2Wdof = CalculateChiSquared(func,ws) + + print 'Chi squared is %s' % chi2 + print 'Chi squared / DOF is %s' % chi2dof + print 'Chi squared weighted is %s' % chi2W + print 'Chi squared weighted / DOF is %s' % chi2Wdof + +Output: + +.. testoutput:: CalculateChiSquaredExample + + Chi squared is 0.0351851851852 + Chi squared / DOF is 0.00439814814815 + Chi squared weighted is 0.0266028783977 + Chi squared weighted / DOF is 0.00332535979971 + + Chi squared is 0.0 + Chi squared / DOF is 0.0 + Chi squared weighted is 0.0 + Chi squared weighted / DOF is 0.0 + +.. categories:: +