Commit 5dbc09a1 authored by Roman Tolchenov's avatar Roman Tolchenov
Browse files

Re #4158. Added FitMW algorithm.

parent e8c2d034
......@@ -83,6 +83,10 @@ public:
std::string parameterDescription(size_t i)const;
/// Checks if a parameter has been set explicitly
bool isExplicitlySet(size_t i)const;
/// Get the fitting error for a parameter
virtual double getError(size_t i) const;
/// Set the fitting error for a parameter
virtual void setError(size_t i, double err);
/// Check if a parameter is active
bool isFixed(size_t i)const;
......@@ -103,6 +107,11 @@ public:
std::string descriptionOfActive(size_t i)const;
/// Check if an active parameter i is actually active
bool isActive(size_t i)const;
/// Return the transformation matrix T between parameters such that p_i = T_ij * ap_j,
/// where ap_j is j-th active parameter.
virtual void getTransformationMatrix(Kernel::Matrix<double>& tm);
/// Is the transformation an identity?
virtual bool isTransformationIdentity() const;
/// Return parameter index from a parameter reference.
size_t getParameterIndex(const ParameterReference& ref)const;
......
......@@ -44,7 +44,9 @@ namespace API
class MANTID_API_DLL FunctionDomain1D: public FunctionDomain
{
public:
FunctionDomain1D(const double startX, const double endX, const size_t n);
FunctionDomain1D(const std::vector<double>& xvalues);
FunctionDomain1D(std::vector<double>::const_iterator from, std::vector<double>::const_iterator to);
/// Return the number of arguments in the domain
virtual size_t size() const {return m_X.size();}
/// Get an x value.
......
......@@ -63,6 +63,8 @@ public:
double* getPointerToCalculated(size_t i);
/// Add other calculated values
FunctionValues& operator+=(const FunctionValues& values);
/// Set all calculated values to zero
void zeroCalculated();
/// set a fitting data value
void setFitData(size_t i,double value);
......
......@@ -8,6 +8,7 @@
#include "MantidAPI/FunctionDomain.h"
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/Jacobian.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/Exception.h"
#include <boost/shared_ptr.hpp>
......@@ -275,6 +276,10 @@ public:
virtual std::string parameterDescription(size_t i)const = 0;
/// Checks if a parameter has been set explicitly
virtual bool isExplicitlySet(size_t i)const = 0;
/// Get the fitting error for a parameter
virtual double getError(size_t i) const = 0;
/// Set the fitting error for a parameter
virtual void setError(size_t i, double err) = 0;
/// Check if a declared parameter i is fixed
virtual bool isFixed(size_t i)const = 0;
......@@ -299,6 +304,11 @@ public:
virtual std::string descriptionOfActive(size_t i)const;
/// Check if an active parameter i is actually active
virtual bool isActive(size_t i)const {return !isFixed(i);}
/// Return the transformation matrix T between parameters such that p_i = T_ij * ap_j,
/// where ap_j is j-th active parameter.
virtual void getTransformationMatrix(Kernel::Matrix<double>& tm);
/// Is the transformation an identity?
virtual bool isTransformationIdentity() const {return true;}
//@}
......@@ -366,6 +376,8 @@ protected:
/// Calculate numerical derivatives
void calNumericalDeriv(const FunctionDomain& domain, Jacobian& out);
/// Calculate the transformation matrix T by numeric differentiation
void calTransformationMatrixNumerically(Kernel::Matrix<double>& tm);
/// Create an instance of a tie without actually tying it to anything
//virtual ParameterTie* createTie(const std::string& parName);
......
......@@ -117,13 +117,13 @@ public:
virtual void function(const FunctionDomain& domain,FunctionValues& values)const;
void functionDeriv(const FunctionDomain& domain, Jacobian& jacobian);
protected:
/// Function you want to fit to.
virtual void function1D(double* out, const double* xValues, const size_t nData)const = 0;
/// Derivatives of function with respect to active parameters
virtual void functionDeriv1D(Jacobian* out, const double* xValues, const size_t nData);
protected:
/// Static reference to the logger class
static Kernel::Logger& g_log;
......@@ -132,6 +132,8 @@ protected:
};
typedef boost::shared_ptr<IFunction1D> IFunction1D_sptr;
} // namespace API
} // namespace Mantid
......
......@@ -88,6 +88,10 @@ public:
virtual std::string parameterDescription(size_t i)const;
/// Checks if a parameter has been set explicitly
virtual bool isExplicitlySet(size_t i)const;
/// Get the fitting error for a parameter
virtual double getError(size_t i) const;
/// Set the fitting error for a parameter
virtual void setError(size_t i, double err);
/// Check if a declared parameter i is active
virtual bool isFixed(size_t i)const;
......@@ -149,6 +153,8 @@ private:
std::vector<std::string> m_parameterNames;
/// Keeps parameter values
std::vector<double> m_parameters;
/// Keeps parameter errors
std::vector<double> m_errors;
/// Holds parameter ties as <parameter index,tie pointer>
std::vector<ParameterTie*> m_ties;
/// Holds the constraints added to function
......
......@@ -123,6 +123,7 @@ std::string CompositeFunction::asString()const
void CompositeFunction::function(const FunctionDomain& domain, FunctionValues& values)const
{
FunctionValues tmp(domain);
values.zeroCalculated();
for(size_t iFun = 0; iFun < nFunctions(); ++iFun)
{
domain.reset();
......@@ -258,6 +259,28 @@ std::string CompositeFunction::parameterDescription(size_t i)const
return ostr.str();
}
/**
* Get the fitting error for a parameter
* @param i :: The index of a parameter
* @return :: the error
*/
double CompositeFunction::getError(size_t i) const
{
size_t iFun = functionIndex(i);
return m_functions[ iFun ]->getError(i - m_paramOffsets[iFun]);
}
/**
* Set the fitting error for a parameter
* @param i :: The index of a parameter
* @param err :: The error value to set
*/
void CompositeFunction::setError(size_t i, double err)
{
size_t iFun = functionIndex(i);
m_functions[ iFun ]->setError(i - m_paramOffsets[iFun], err);
}
/// Value of i-th active parameter. Override this method to make fitted parameters different from the declared
double CompositeFunction::activeParameter(size_t i)const
{
......@@ -837,5 +860,50 @@ IFunction_sptr CompositeFunction::getContainingFunction(const ParameterReference
// return IFunction_sptr();
//}
/// Is the transformation an identity?
bool CompositeFunction::isTransformationIdentity() const
{
for(size_t iFun = 0; iFun < nFunctions(); ++iFun)
{
if ( !getFunction(iFun)->isTransformationIdentity() )
{
return false;
}
}
return true;
}
/**
* Return the transformation matrix T between parameters such that p_i = T_ij * ap_j,
* where ap_j is j-th active parameter.
* @param tm :: The output matrix
*/
void CompositeFunction::getTransformationMatrix(Kernel::Matrix<double>& tm)
{
size_t np = nParams();
tm.setMem(np, np);
tm.identityMatrix();
if ( this->isTransformationIdentity() ) return;
for(size_t iFun = 0; iFun < nFunctions(); ++iFun)
{
IFunction_sptr f = getFunction(iFun);
if ( !f->isTransformationIdentity() )
{
size_t i0 = m_paramOffsets[iFun];
size_t np1 = f->nParams();
Kernel::Matrix<double> tm1(np1,np1);
f->getTransformationMatrix(tm1);
for(size_t i = 0; i < np1; ++i)
{
for(size_t j = 0; j < np1; ++j)
{
tm[i0 + i][i0 + j] = tm1[i][j];
}
}
} // !identity
}
}
} // namespace API
} // namespace Mantid
......@@ -22,6 +22,35 @@ FunctionDomain1D::FunctionDomain1D(const std::vector<double>& xvalues)
m_X.assign(xvalues.begin(),xvalues.end());
}
FunctionDomain1D::FunctionDomain1D(std::vector<double>::const_iterator from, std::vector<double>::const_iterator to)
{
if (from == to)
{
throw std::invalid_argument("FunctionDomain1D cannot have zero size.");
}
m_X.assign(from,to);
}
FunctionDomain1D::FunctionDomain1D(const double startX, const double endX, const size_t n)
{
if (n == 0)
{
throw std::invalid_argument("FunctionDomain1D cannot have zero size.");
}
m_X.resize(n);
if (n == 1)
{
m_X[0] = (startX + endX) / 2;
}
else
{
const double dx = (endX - startX) / (n - 1);
for(size_t i = 0; i < n; ++i)
{
m_X[i] = startX + dx * i;
}
}
}
} // namespace API
} // namespace Mantid
......@@ -66,6 +66,13 @@ namespace API
return *this;
}
/// Set all calculated values to zero
void FunctionValues::zeroCalculated()
{
std::fill(m_calculated.begin(),m_calculated.end(),0.0);
}
/**
* Set a fitting data value.
* @param i :: Index
......
......@@ -12,6 +12,7 @@
#include <Poco/StringTokenizer.h>
#include <limits>
#include <sstream>
#include <iostream>
......@@ -457,7 +458,7 @@ std::string IFunction::descriptionOfActive(size_t i)const
void IFunction::calNumericalDeriv(const FunctionDomain& domain, Jacobian& jacobian)
{
const double minDouble = std::numeric_limits<double>::min();
const double epsilon = std::numeric_limits<double>::epsilon();
const double epsilon = std::numeric_limits<double>::epsilon() * 100;
double stepPercentage = 0.001; // step percentage
double step; // real step
double cutoff = 100.0*minDouble/stepPercentage;
......@@ -502,5 +503,66 @@ void IFunction::calNumericalDeriv(const FunctionDomain& domain, Jacobian& jacobi
}
}
/**
* Return the transformation matrix T between parameters such that p_i = T_ij * ap_j,
* where ap_j is j-th active parameter.
* @param tm :: The output matrix
*/
void IFunction::getTransformationMatrix(Kernel::Matrix<double>& tm)
{
size_t np = nParams();
tm.setMem(np, np);
tm.identityMatrix();
}
/**
* Calculate the transformation matrix T by numeric differentiation
* @param tm :: The output transformation matrix.
*/
void IFunction::calTransformationMatrixNumerically(Kernel::Matrix<double>& tm)
{
double epsilon = std::numeric_limits<double>::epsilon() * 100;
size_t np = nParams();
tm.setMem(np,np);
for(size_t i = 0; i < np; ++i)
{
double p0 = getParameter(i);
for(size_t j = 0; j < np; ++j)
{
double ap = activeParameter(j);
double step = ap == 0.0? epsilon : ap * epsilon;
setActiveParameter( j, ap + step );
tm[i][j] = ( getParameter(i) - p0 ) / step;
setActiveParameter( j, ap );
}
}
}
} // namespace API
} // namespace Mantid
///\cond TEMPLATE
namespace Mantid
{
namespace Kernel
{
template<> MANTID_API_DLL
boost::shared_ptr<Mantid::API::IFunction> IPropertyManager::getValue<boost::shared_ptr<Mantid::API::IFunction> >(const std::string &name) const
{
PropertyWithValue<boost::shared_ptr<Mantid::API::IFunction> >* prop =
dynamic_cast<PropertyWithValue<boost::shared_ptr<Mantid::API::IFunction> >*>(getPointerToProperty(name));
if (prop)
{
return *prop;
}
else
{
std::string message = "Attempt to assign property "+ name +" to incorrect type. Expected IFitFunction.";
throw std::runtime_error(message);
}
}
} // namespace Kernel
} // namespace Mantid
///\endcond TEMPLATE
......@@ -204,6 +204,35 @@ std::string ParamFunction::parameterDescription(size_t i)const
return m_parameterDescriptions[i];
}
/**
* Get the fitting error for a parameter
* @param i :: The index of a parameter
* @return :: the error
*/
double ParamFunction::getError(size_t i) const
{
if (i >= nParams())
{
throw std::out_of_range("ParamFunction parameter index out of range.");
}
return m_errors[i];
}
/**
* Set the fitting error for a parameter
* @param i :: The index of a parameter
* @param err :: The error value to set
*/
void ParamFunction::setError(size_t i, double err)
{
if (i >= nParams())
{
throw std::out_of_range("ParamFunction parameter index out of range.");
}
m_errors[i] = err;
}
/**
* Declare a new parameter. To used in the implementation'c constructor.
* @param name :: The parameter name.
......@@ -227,6 +256,7 @@ void ParamFunction::declareParameter(const std::string& name,double initValue, c
m_parameterNames.push_back(ucName);
m_parameterDescriptions.push_back(description);
m_parameters.push_back(initValue);
m_errors.push_back(0.0);
m_explicitlySet.push_back(false);
}
......
......@@ -4,9 +4,9 @@ set ( SRC_FILES
src/BackToBackExponential.cpp
src/BackgroundFunction.cpp
# src/BivariateNormal.cpp
# src/BoundaryConstraint.cpp
src/BoundaryConstraint.cpp
src/Chebyshev.cpp
# src/Convolution.cpp
src/Convolution.cpp
src/CostFuncFitting.cpp
# src/CostFuncIgnorePosPeaks.cpp
src/CostFuncLeastSquares.cpp
......@@ -17,28 +17,28 @@ set ( SRC_FILES
src/ExpDecayMuon.cpp
src/ExpDecayOsc.cpp
src/FlatBackground.cpp
# src/FRConjugateGradientMinimizer.cpp
src/FRConjugateGradientMinimizer.cpp
# src/Fit.cpp
# src/Fit1D.cpp
src/Fit1D.cpp
src/FitMW.cpp
src/FuncMinimizerFactory.cpp
src/GausDecay.cpp
src/GausOsc.cpp
src/Gaussian.cpp
src/MuonFInteraction.cpp
# src/Gaussian1D.cpp
# src/GaussianLinearBG1D.cpp
src/Gaussian1D.cpp
src/GaussianLinearBG1D.cpp
# src/GenericFit.cpp
src/IFuncMinimizer.cpp
# src/NewFit.cpp
# src/IkedaCarpenterPV.cpp
# src/LevenbergMarquardtMinimizer.cpp
src/LevenbergMarquardtMinimizer.cpp
# src/Linear.cpp
src/LinearBackground.cpp
src/LogNormal.cpp
src/Lorentzian.cpp
# src/Lorentzian1D.cpp
src/Lorentzian1D.cpp
# src/MultiBG.cpp
# src/PRConjugateGradientMinimizer.cpp
src/PRConjugateGradientMinimizer.cpp
# src/PlotPeakByLogValue.cpp
# src/ProductFunctionMW.cpp
src/Quadratic.cpp
......@@ -47,12 +47,13 @@ set ( SRC_FILES
src/SCDPanelErrors.cpp
src/StaticKuboToyabe.cpp
src/SimplexMinimizer.cpp
src/SteepestDescentMinimizer.cpp
# src/SpecialFunctionHelper.cpp
src/SplineBackground.cpp
src/StretchExp.cpp
src/StretchExpMuon.cpp
src/UserFunction.cpp
# src/UserFunction1D.cpp
src/UserFunction1D.cpp
)
set ( SRC_UNITY_IGNORE_FILES src/Fit1D.cpp src/GSLFunctions.cpp )
......@@ -63,9 +64,9 @@ set ( INC_FILES
inc/MantidCurveFitting/BackToBackExponential.h
inc/MantidCurveFitting/BackgroundFunction.h
# inc/MantidCurveFitting/BivariateNormal.h
# inc/MantidCurveFitting/BoundaryConstraint.h
inc/MantidCurveFitting/BoundaryConstraint.h
inc/MantidCurveFitting/Chebyshev.h
# inc/MantidCurveFitting/Convolution.h
inc/MantidCurveFitting/Convolution.h
inc/MantidCurveFitting/CostFuncFitting.h
# inc/MantidCurveFitting/CostFuncIgnorePosPeaks.h
inc/MantidCurveFitting/CostFuncLeastSquares.h
......@@ -77,22 +78,24 @@ set ( INC_FILES
inc/MantidCurveFitting/ExpDecayMuon.h
inc/MantidCurveFitting/ExpDecayOsc.h
inc/MantidCurveFitting/FlatBackground.h
# inc/MantidCurveFitting/FRConjugateGradientMinimizer.h
inc/MantidCurveFitting/FRConjugateGradientMinimizer.h
# inc/MantidCurveFitting/Fit.h
# inc/MantidCurveFitting/Fit1D.h
inc/MantidCurveFitting/Fit1D.h
inc/MantidCurveFitting/FitMW.h
inc/MantidCurveFitting/FuncMinimizerFactory.h
inc/MantidCurveFitting/Jacobian.h
inc/MantidCurveFitting/GSLJacobian.h
inc/MantidCurveFitting/GSLMatrix.h
inc/MantidCurveFitting/GSLVector.h
inc/MantidCurveFitting/GausDecay.h
inc/MantidCurveFitting/GausOsc.h
inc/MantidCurveFitting/Gaussian.h
# inc/MantidCurveFitting/Gaussian1D.h
# inc/MantidCurveFitting/GaussianLinearBG1D.h
# inc/MantidCurveFitting/GenericFit.h
# inc/MantidCurveFitting/NewFit.h
inc/MantidCurveFitting/Gaussian1D.h
inc/MantidCurveFitting/GaussianLinearBG1D.h
inc/MantidCurveFitting/GenericFit.h
inc/MantidCurveFitting/IFuncMinimizer.h
inc/MantidCurveFitting/IkedaCarpenterPV.h
# inc/MantidCurveFitting/LevenbergMarquardtMinimizer.h
inc/MantidCurveFitting/LevenbergMarquardtMinimizer.h
# inc/MantidCurveFitting/Linear.h
inc/MantidCurveFitting/LinearBackground.h
inc/MantidCurveFitting/LogNormal.h
......@@ -100,7 +103,7 @@ set ( INC_FILES
inc/MantidCurveFitting/Lorentzian1D.h
inc/MantidCurveFitting/MuonFInteraction.h
# inc/MantidCurveFitting/MultiBG.h
# inc/MantidCurveFitting/PRConjugateGradientMinimizer.h
inc/MantidCurveFitting/PRConjugateGradientMinimizer.h
# inc/MantidCurveFitting/PlotPeakByLogValue.h
# inc/MantidCurveFitting/ProductFunctionMW.h
inc/MantidCurveFitting/Quadratic.h
......@@ -109,12 +112,13 @@ set ( INC_FILES
inc/MantidCurveFitting/SCDPanelErrors.h
inc/MantidCurveFitting/StaticKuboToyabe.h
inc/MantidCurveFitting/SimplexMinimizer.h
inc/MantidCurveFitting/SteepestDescentMinimizer.h
# inc/MantidCurveFitting/SpecialFunctionSupport.h
inc/MantidCurveFitting/SplineBackground.h
inc/MantidCurveFitting/StretchExp.h
inc/MantidCurveFitting/StretchExpMuon.h
inc/MantidCurveFitting/UserFunction.h
# inc/MantidCurveFitting/UserFunction1D.h
inc/MantidCurveFitting/UserFunction1D.h
)
set ( TEST_FILES
......@@ -130,30 +134,33 @@ set ( TEST_FILES
# test/BackToBackExponentialTest.h
test/BFGSTest.h
# test/BivariateNormalTest.h
# test/BoundaryConstraintTest.h
test/BoundaryConstraintTest.h
# test/ChebyshevTest.h
# test/CompositeFunctionTest.h
test/CompositeFunctionTest.h
# test/ConvolutionTest.h
# test/ExpDecayTest.h
test/ExpDecayTest.h
# test/ExpDecayMuonTest.h
# test/ExpDecayOscTest.h
# test/FitTest.h
test/FitMWTest.h
test/FRConjugateGradientTest.h
# test/NewFitTest.h
# test/FuncMinimizerFactoryTest.h
# test/FunctionTest.h
# test/FunctionFactoryTest.h
# test/Gaussian1DTest.h
# test/GausDecayTest.h
# test/GaussianTest.h
test/GausDecayTest.h
test/GaussianTest.h
# test/IkedaCarpenterPVTest.h
test/LeastSquaresTest.h
# test/LinearBackgroundTest.h
test/LinearBackgroundTest.h
# test/LinearTest.h
# test/LogNormalTest.h
# test/Lorentzian1DTest.h
# test/LorentzianTest.h
test/LorentzianTest.h
# test/MultiBGTest.h
# test/PlotPeakByLogValueTest.h
test/PRConjugateGradientTest.h
# test/ProductFunctionMWTest.h
# test/QuadraticBackgroundTest.h
# test/QuadraticTest.h
......@@ -186,9 +193,9 @@ include_directories ( inc )
target_link_libraries ( CurveFitting ${MANTIDLIBS} ${GSL_LIBRARIES} )
if ( CXXTEST_FOUND )
include_directories ( ../DataHandling/inc )
include_directories ( ../DataHandling/inc ../TestHelpers/inc )
cxxtest_add_test ( CurveFittingTest ${TEST_FILES} )
target_link_libraries( CurveFittingTest CurveFitting DataHandling )
target_link_libraries( CurveFittingTest CurveFitting DataHandling DataObjects)
add_dependencies ( FrameworkTests CurveFittingTest )
# Add to the 'FrameworkTests' group in VS
set_property ( TARGET CurveFittingTest PROPERTY FOLDER "UnitTests" )
......
......@@ -123,11 +123,17 @@ namespace Mantid
virtual const std::string category() const { return "General";}
/// Writes itself into a string
std::string asString()const;
void functionMW(double* out, const double* xValues, const size_t nData)const;
void functionDerivMW(API::Jacobian* out, const double* xValues, const size_t nData);
/// Function you want to fit to.
/// @param domain :: The buffer for writing the calculated values. Must be big enough to accept dataSize() values
virtual void function(const API::FunctionDomain& domain, API::FunctionValues& values)const;
/// Derivatives of function with respect to active parameters
virtual void functionDeriv(const API::FunctionDomain& domain, API::Jacobian& jacobian);
void function1D(double* out, const double* xValues, const size_t nData)const;
void functionDeriv1D(API::Jacobian* out, const double* xValues, const size_t nData);
/// Add a function.
size_t addFunction(IFunction* f);
size_t addFunction(API::IFunction_sptr f);
/// Deletes and zeroes pointer m_resolution forsing function(...) to recalculate the resolution function
void refreshResolution()const;
......
......@@ -6,8 +6,9 @@
//----------------------------------------------------------------------
#include "MantidAPI/ICostFunction.h"
#include "MantidAPI/IFunction.h"
#include "MantidKernel/Matrix.h"
#include <gsl/gsl_matrix.h>
//#include <gsl/gsl_matrix.h>