Commit ca899c5e authored by Roman Tolchenov's avatar Roman Tolchenov
Browse files

Re #15700. Function almost ready for fitting.

parent d4eb23d3
......@@ -149,6 +149,8 @@ public:
bool removeTie(size_t i) override;
/// Get the tie of i-th parameter
ParameterTie *getTie(size_t i) const override;
/// Add a new tie
void addTie(ParameterTie *tie) override;
/// Overwrite IFunction methods
void addConstraint(IConstraint *ic) override;
......@@ -181,6 +183,8 @@ public:
std::string parameterLocalName(size_t i) const;
/// Check the function.
void checkFunction();
/// Remove all member functions
void clear();
/// Returns the number of attributes associated with the function
virtual size_t nLocalAttributes() const { return 0; }
......@@ -220,8 +224,6 @@ protected:
/// Declare a new parameter
void declareParameter(const std::string &name, double initValue = 0,
const std::string &description = "") override;
/// Add a new tie
void addTie(ParameterTie *tie) override;
size_t paramOffset(size_t i) const { return m_paramOffsets[i]; }
......@@ -243,7 +245,6 @@ private:
size_t m_nParams;
/// Function counter to be used in nextConstraint
mutable size_t m_iConstraintFunction;
/// Flag set to use numerical derivatives
};
/// shared pointer to the composite function base class
......
......@@ -464,6 +464,8 @@ public:
virtual bool removeTie(size_t i) = 0;
/// Get the tie of i-th parameter
virtual ParameterTie *getTie(size_t i) const = 0;
/// Add a new tie. Derived classes must provide storage for ties
virtual void addTie(ParameterTie *tie) = 0;
//@}
/** @name Constraints */
......@@ -543,9 +545,6 @@ protected:
boost::shared_ptr<const MatrixWorkspace> ws,
size_t wsIndex) const;
/// Add a new tie. Derived classes must provide storage for ties
virtual void addTie(ParameterTie *tie) = 0;
/// Override to declare function attributes
virtual void declareAttributes() {}
/// Override to declare function parameters
......
......@@ -110,6 +110,8 @@ public:
bool removeTie(size_t i) override;
/// Get the tie of i-th parameter
ParameterTie *getTie(size_t i) const override;
/// Add a new tie
void addTie(ParameterTie *tie) override;
/// Add a constraint to function
void addConstraint(IConstraint *ic) override;
......@@ -125,8 +127,6 @@ protected:
void declareParameter(const std::string &name, double initValue = 0,
const std::string &description = "") override;
/// Add a new tie
void addTie(ParameterTie *tie) override;
/// Get the address of the parameter. For use in UserFunction with mu::Parser
virtual double *getParameterAddress(size_t i);
......
......@@ -375,6 +375,16 @@ void CompositeFunction::checkFunction() {
}
}
/**
* Remove all member functions
*/
void CompositeFunction::clear() {
m_nParams = 0;
m_paramOffsets.clear();
m_IFunction.clear();
m_functions.clear();
}
/** Add a function
* @param f :: A pointer to the added function
* @return The function index
......@@ -736,5 +746,6 @@ CompositeFunction::getContainingFunction(const ParameterReference &ref) const {
return IFunction_sptr();
}
} // namespace API
} // namespace Mantid
......@@ -5,7 +5,7 @@
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/CompositeFunction.h"
#include "MantidAPI/IFunction1D.h"
#include "MantidAPI/IFunction.h"
#include "MantidCurveFitting/Functions/CrystalFieldPeaks.h"
namespace Mantid {
......@@ -35,7 +35,7 @@ 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 CrystalFieldSpectrum : public API::IFunction1D {
class DLLExport CrystalFieldSpectrum : public API::IFunction {
public:
/// Constructor
......@@ -90,6 +90,9 @@ public:
/// Return parameter index from a parameter reference.
virtual size_t getParameterIndex(const API::ParameterReference &ref) const override;
/// Tie a parameter to other parameters (or a constant)
virtual API::ParameterTie *tie(const std::string &parName, const std::string &expr,
bool isDefault = false) override;
/// Apply the ties
virtual void applyTies() override;
/// Remove all ties
......@@ -109,6 +112,9 @@ public:
/// Set up the function for a fit.
virtual void setUpForFit() override;
/// Build m_spectrum function.
void buildSpectrumFunction() const;
protected:
/// Declare a new parameter
virtual void declareParameter(const std::string &name, double initValue = 0,
......@@ -118,25 +124,41 @@ protected:
virtual void addTie(API::ParameterTie *tie) override;
//@}
/// Evaluate the function
void function1D(double *out, const double *xValues,
const size_t nData) const override;
public:
/** @name Attributes */
//@{
/// Returns the number of attributes associated with the function
virtual size_t nAttributes() const override;
/// Returns a list of attribute names
virtual std::vector<std::string> getAttributeNames() const override;
/// Return a value of attribute attName
virtual Attribute getAttribute(const std::string &name) const override;
/// Set a value to attribute attName
virtual void setAttribute(const std::string &name, const Attribute &) override;
/// Check if attribute attName exists
virtual bool hasAttribute(const std::string &name) const override;
//@}
void setAttribute(const std::string &attName, const Attribute &) override;
/// Evaluate the function
void function(const API::FunctionDomain &domain, API::FunctionValues &values) const override;
protected:
/// overwrite IFunction base class method, which declare function parameters
void init() override;
private:
/// Build m_spectrum function.
void buildSpectrumFunction();
/// Test if a name (parameter's or attribute's) belongs to m_crystalFiled
bool isOwnName(const std::string &aName) const;
/// Update m_spectrum function.
void updateSpectrumFunction() const;
/// Function that calculates peak centres and intensities.
CrystalFieldPeaks m_crystalField;
/// Function that actually callculates the spectrum.
API::CompositeFunction m_spectrum;
mutable API::CompositeFunction m_spectrum;
/// Cached number of parameters in m_crystalField.
size_t m_nOwnParams;
/// Flag indicating that updateSpectrumFunction() is required.
mutable bool m_dirty;
};
} // namespace Functions
......
......@@ -50,9 +50,15 @@ public:
double centre() const override { return getParameter("PeakCentre"); }
double height() const override;
double fwhm() const override { return getParameter("FWHM"); }
double intensity() const override { return getParameter("Amplitude"); }
void setCentre(const double c) override { setParameter("PeakCentre", c); }
void setHeight(const double h) override;
void setFwhm(const double w) override { setParameter("FWHM", w); }
void setIntensity(const double i) override { setParameter("Amplitude", i); }
void fixCentre() override;
void unfixCentre() override;
void fixIntensity() override;
void unfixIntensity() override;
/// overwrite IFunction base class methods
std::string name() const override { return "Lorentzian"; }
......
......@@ -3,8 +3,9 @@
//----------------------------------------------------------------------
#include "MantidCurveFitting/Functions/CrystalFieldSpectrum.h"
#include "MantidCurveFitting/Functions/DeltaFunction.h"
#include "MantidAPI/IFunction1D.h"
#include "MantidAPI/IFunction.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/ParameterTie.h"
#include <cmath>
#include <algorithm>
......@@ -31,9 +32,9 @@ using namespace API;
DECLARE_FUNCTION(CrystalFieldSpectrum)
/// Constructor
CrystalFieldSpectrum::CrystalFieldSpectrum(): m_nOwnParams(m_crystalField.nParams()) {
CrystalFieldSpectrum::CrystalFieldSpectrum(): m_nOwnParams(m_crystalField.nParams()),m_dirty(true) {
declareAttribute("PeakShape", Attribute("Lorentzian"));
setAttributeValue("NumDeriv", true);
declareAttribute("FWHM", Attribute(0.0));
}
void CrystalFieldSpectrum::init() {}
......@@ -43,7 +44,11 @@ void CrystalFieldSpectrum::setParameter(size_t i, const double &value,
bool explicitlySet) {
if (i < m_nOwnParams) {
m_crystalField.setParameter(i, value, explicitlySet);
m_dirty = true;
} else {
if (m_dirty) {
updateSpectrumFunction();
}
m_spectrum.setParameter(i - m_nOwnParams, value, explicitlySet);
}
}
......@@ -54,6 +59,9 @@ void CrystalFieldSpectrum::setParameterDescription(size_t i,
if (i < m_nOwnParams) {
m_crystalField.setParameterDescription(i, description);
} else {
if (m_dirty) {
updateSpectrumFunction();
}
m_spectrum.setParameterDescription(i - m_nOwnParams, description);
}
}
......@@ -91,13 +99,13 @@ size_t CrystalFieldSpectrum::nParams() const {
/// Returns the index of parameter name
size_t CrystalFieldSpectrum::parameterIndex(const std::string &name) const {
if (name.empty()) {
throw std::invalid_argument("Parameter name cannot be empty string.");
}
if (name.front() != 'f' || name.find('.') == std::string::npos) {
if (isOwnName(name)) {
return m_crystalField.parameterIndex(name);
} else {
return m_spectrum.parameterIndex(name);
if (m_dirty) {
updateSpectrumFunction();
}
return m_spectrum.parameterIndex(name) + m_nOwnParams;
}
}
......@@ -126,6 +134,9 @@ void CrystalFieldSpectrum::setError(size_t i, double err) {
if (i < m_nOwnParams) {
m_crystalField.setError(i, err);
} else {
if (m_dirty) {
updateSpectrumFunction();
}
m_spectrum.setError(i - m_nOwnParams, err);
}
}
......@@ -158,7 +169,20 @@ size_t CrystalFieldSpectrum::getParameterIndex(const ParameterReference &ref) co
if (index < m_nOwnParams) {
return index;
}
return m_spectrum.getParameterIndex(ref);
return m_spectrum.getParameterIndex(ref) + m_nOwnParams;
}
/// Tie a parameter to other parameters (or a constant)
API::ParameterTie *CrystalFieldSpectrum::tie(const std::string &parName, const std::string &expr,
bool isDefault) {
if (isOwnName(parName)) {
return m_crystalField.tie(parName, expr, isDefault);
} else {
if (m_dirty) {
updateSpectrumFunction();
}
return m_spectrum.tie(parName, expr, isDefault);
}
}
/// Apply the ties
......@@ -210,6 +234,7 @@ void CrystalFieldSpectrum::removeConstraint(const std::string &parName) {
/// Set up the function for a fit.
void CrystalFieldSpectrum::setUpForFit() {
updateSpectrumFunction();
}
/// Declare a new parameter
......@@ -220,29 +245,164 @@ void CrystalFieldSpectrum::declareParameter(const std::string &, double,
}
/// Add a new tie. Derived classes must provide storage for ties
void CrystalFieldSpectrum::addTie(ParameterTie *tie) {
void CrystalFieldSpectrum::addTie(API::ParameterTie *tie) {
size_t i = getParameterIndex(*tie);
if (i < m_nOwnParams) {
m_crystalField.addTie(tie);
} else {
if (m_dirty) {
updateSpectrumFunction();
}
m_spectrum.addTie(tie);
}
}
/// Returns the number of attributes associated with the function
size_t CrystalFieldSpectrum::nAttributes() const {
return IFunction::nAttributes() + m_crystalField.nAttributes() + m_spectrum.nAttributes();
}
/// Returns a list of attribute names
std::vector<std::string> CrystalFieldSpectrum::getAttributeNames() const {
std::vector<std::string> attNames(2);
attNames[0] = "PeakShape";
attNames[1] = "FWHM";
auto cfNames = m_crystalField.getAttributeNames();
auto spNames = m_spectrum.getAttributeNames();
attNames.insert(attNames.end(), cfNames.begin(), cfNames.end());
attNames.insert(attNames.end(), spNames.begin(), spNames.end());
return attNames;
}
/// Return a value of attribute attName
API::IFunction::Attribute CrystalFieldSpectrum::getAttribute(const std::string &attName) const {
if (attName == "PeakShape" || attName == "FWHM") {
return IFunction::getAttribute(attName);
} else if (attName == "NumDeriv") {
return m_spectrum.getAttribute(attName);
} else if (isOwnName(attName)) {
return m_crystalField.getAttribute(attName);
} else {
return m_spectrum.getAttribute(attName);
}
}
/// Set a value to attribute attName
void CrystalFieldSpectrum::setAttribute(const std::string &attName,
const IFunction::Attribute &att) {
if (attName == "PeakShape") {
if (attName == "PeakShape" || attName == "FWHM") {
IFunction::setAttribute(attName, att);
m_dirty = true;
m_spectrum.clear();
} else if (attName == "NumDeriv") {
m_spectrum.setAttribute(attName, att);
} else if (isOwnName(attName)) {
m_crystalField.setAttribute(attName, att);
m_dirty = true;
} else {
if (m_dirty) {
updateSpectrumFunction();
}
m_spectrum.setAttribute(attName, att);
}
}
/// Check if attribute attName exists
bool CrystalFieldSpectrum::hasAttribute(const std::string &attName) const {
if (attName == "NumDeriv" || IFunction::hasAttribute(attName)) {
return true;
}
if (isOwnName(attName)) {
return m_crystalField.hasAttribute(attName);
} else {
return m_spectrum.hasAttribute(attName);
}
//CompositeFunction::setAttribute(attName, att);
}
// Evaluates the function
void CrystalFieldSpectrum::function1D(double *out, const double *xValues,
const size_t nData) const {
void CrystalFieldSpectrum::function(const API::FunctionDomain &domain,
API::FunctionValues &values) const {
updateSpectrumFunction();
m_spectrum.function(domain, values);
}
/// Test if a name (parameter's or attribute's) belongs to m_crystalFiled
/// @param aName :: A name to test.
bool CrystalFieldSpectrum::isOwnName(const std::string &aName) const {
if (aName.empty()) {
throw std::invalid_argument("Parameter or attribute name cannot be empty string.");
}
return (aName.front() != 'f' || aName.find('.') == std::string::npos);
}
/// Uses m_crystalField to calculate peak centres and intensities
/// then populates m_spectrum with peaks of type given in PeakShape attribute.
void CrystalFieldSpectrum::buildSpectrumFunction() {
void CrystalFieldSpectrum::buildSpectrumFunction() const {
m_dirty = false;
m_spectrum.clear();
FunctionDomainGeneral domain;
FunctionValues values;
m_crystalField.functionGeneral(domain, values);
if (values.size() == 0) {
return;
}
if (values.size() % 2 != 0) {
throw std::runtime_error("CrystalFieldPeaks returned add number of values.");
}
auto peakShape = IFunction::getAttribute("PeakShape").asString();
auto fwhm = IFunction::getAttribute("FWHM").asDouble();
auto nPeaks = values.size() / 2;
for(size_t i = 0; i < nPeaks; ++i) {
auto fun = API::FunctionFactory::Instance().createFunction(peakShape);
auto peak = boost::dynamic_pointer_cast<API::IPeakFunction>(fun);
if (!peak) {
throw std::runtime_error("A peak function is expected.");
}
peak->fixCentre();
peak->fixIntensity();
peak->setCentre(values.getCalculated(i));
peak->setIntensity(values.getCalculated(i + nPeaks));
peak->setIntensity(fwhm);
m_spectrum.addFunction(peak);
}
}
/// Update m_spectrum function.
void CrystalFieldSpectrum::updateSpectrumFunction() const {
if (m_spectrum.nFunctions() == 0) {
buildSpectrumFunction();
return;
}
m_dirty = false;
FunctionDomainGeneral domain;
FunctionValues values;
m_crystalField.functionGeneral(domain, values);
auto nPeaks = values.size() / 2;
if (m_spectrum.nFunctions() != nPeaks) {
//throw std::runtime_error("Number of peaks has changed: "+std::to_string(nPeaks));
buildSpectrumFunction();
return;
}
for(size_t i = 0; i < nPeaks; ++i) {
auto fun = m_spectrum.getFunction(i);
auto peak = boost::dynamic_pointer_cast<API::IPeakFunction>(fun);
if (!peak) {
throw std::runtime_error("A peak function is expected.");
}
auto centre = values.getCalculated(i);
auto intensity = values.getCalculated(i + nPeaks);
peak->setCentre(centre);
peak->setIntensity(intensity);
}
}
} // namespace Functions
} // namespace CurveFitting
......
......@@ -41,6 +41,22 @@ void Lorentzian::setHeight(const double h) {
}
}
void Lorentzian::fixCentre() {
fixParameter("PeakCentre");
}
void Lorentzian::unfixCentre() {
unfixParameter("PeakCentre");
}
void Lorentzian::fixIntensity() {
fixParameter("Amplitude");
}
void Lorentzian::unfixIntensity() {
unfixParameter("Amplitude");
}
void Lorentzian::functionLocal(double *out, const double *xValues,
const size_t nData) const {
const double amplitude = getParameter("Amplitude");
......
......@@ -3,7 +3,11 @@
#include <cxxtest/TestSuite.h>
#include "MantidAPI/AlgorithmFactory.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidCurveFitting/Functions/CrystalFieldSpectrum.h"
#include "MantidCurveFitting/Algorithms/Fit.h"
......@@ -16,15 +20,112 @@ using namespace Mantid::CurveFitting;
using namespace Mantid::CurveFitting::Algorithms;
using namespace Mantid::CurveFitting::Functions;
//typedef Mantid::DataObjects::Workspace2D_sptr WS_type;
//typedef Mantid::DataObjects::TableWorkspace_sptr TWS_type;
class CrystalFieldSpectrumTest : public CxxTest::TestSuite {
public:
void testFunction() {
void xtest_function() {
std::cerr << "\n\nTest function\n\n";
CrystalFieldSpectrum fun;
fun.setParameter("B20", 0.37737);
fun.setParameter("B22", 3.9770);
fun.setParameter("B40", -0.031787);
fun.setParameter("B42", -0.11611);
fun.setParameter("B44", -0.12544);
fun.setAttributeValue("Ion", "Ce");
fun.setAttributeValue("Temperature", 44.0);
fun.setAttributeValue("ToleranceIntensity", 0.001);
fun.buildSpectrumFunction();
std::cerr << fun.nParams() << ' ' << fun.nAttributes() << std::endl;
auto attNames = fun.getAttributeNames();
auto parNames = fun.getParameterNames();
TS_ASSERT_EQUALS(fun.nAttributes(), attNames.size());
TS_ASSERT_EQUALS(fun.nParams(), parNames.size());
//std::cerr << "\nAttributes:" << std::endl;
//for(auto &name: attNames) {
// std::cerr << name << std::endl;
//}
//std::cerr << "\nParameters:" << std::endl;
//for(size_t i = 0; i < parNames.size(); ++i) {
// auto &name = parNames[i];
// std::cerr << i << ' ' << fun.parameterIndex(name) << ' ' << name << ' ' << fun.isFixed(i) << fun.isActive(i) << std::endl;
//}
auto i = fun.parameterIndex("f0.Amplitude");
TS_ASSERT(fun.isFixed(i));
TS_ASSERT(!fun.isActive(i));
i = fun.parameterIndex("f0.PeakCentre");
TS_ASSERT(fun.isFixed(i));
TS_ASSERT(!fun.isActive(i));
i = fun.parameterIndex("f0.FWHM");
TS_ASSERT(!fun.isFixed(i));
TS_ASSERT(fun.isActive(i));
i = fun.parameterIndex("f1.Amplitude");
TS_ASSERT(fun.isFixed(i));
TS_ASSERT(!fun.isActive(i));
i = fun.parameterIndex("f1.PeakCentre");
TS_ASSERT(fun.isFixed(i));
TS_ASSERT(!fun.isActive(i));
i = fun.parameterIndex("f1.FWHM");
TS_ASSERT(!fun.isFixed(i));
TS_ASSERT(fun.isActive(i));
i = fun.parameterIndex("f2.Amplitude");
TS_ASSERT(fun.isFixed(i));
TS_ASSERT(!fun.isActive(i));
i = fun.parameterIndex("f2.PeakCentre");
TS_ASSERT(fun.isFixed(i));
TS_ASSERT(!fun.isActive(i));
i = fun.parameterIndex("f2.FWHM");
TS_ASSERT(!fun.isFixed(i));
TS_ASSERT(fun.isActive(i));
std::cerr << fun.getParameter("f0.PeakCentre") << ' ' << fun.getParameter("f0.Amplitude") << ' ' << fun.getParameter("f0.FWHM") << std::endl;
std::cerr << fun.getParameter("f1.PeakCentre") << ' ' << fun.getParameter("f1.Amplitude") << ' ' << fun.getParameter("f1.FWHM") << std::endl;
std::cerr << fun.getParameter("f2.PeakCentre") << ' ' << fun.getParameter("f2.Amplitude") << ' ' << fun.getParameter("f2.FWHM") << std::endl;
}
void xtest_evaluate() {
auto fun = makeFunction();
fun->setParameter("f0.FWHM", 2.0);
fun->setParameter("f1.FWHM", 2.0);
fun->setParameter("f2.FWHM", 2.0);
FunctionDomain1DVector x(0.0, 50.0, 100);
FunctionValues y(x);
fun->function(x, y);
for(size_t i = 0; i < x.size(); ++i) {
std::cerr << x[i] << ' ' << y[i] << std::endl;
}
}