Mantid
Loading...
Searching...
No Matches
FickDiffusionSQE.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7// Mantid Coding standars <http://www.mantidproject.org/Coding_Standards>
8// Main Module Header
10// Mantid Headers from the same project
11// N/A
12// Mantid headers from other projects
14#include "MantidAPI/IFunction.h"
15#include "MantidAPI/Jacobian.h"
16// third party library headers
17// N/A
18// standard library headers
19#include <cmath>
20#include <limits>
21
22namespace {
23Mantid::Kernel::Logger g_log("FickDiffusionSQE");
24}
25
27
28DECLARE_FUNCTION(FickDiffusionSQE)
29
30
33void FickDiffusionSQE::declareParameters() {
34 this->declareParameter("Height", 1.0, "scaling factor");
35 this->declareParameter("Centre", 0.0, "Shift along the X-axis");
36 this->declareParameter("DiffCoeff", 2.3, "Diffusion coefficient (10^(-5)cm^2/s)");
37}
38
45void FickDiffusionSQE::function1D(double *out, const double *xValues, const size_t nData) const {
46 auto H = this->getParameter("Height");
47 auto C = this->getParameter("Centre");
48 auto Q = this->getAttribute("Q").asDouble();
49 auto D = this->getParameter("DiffCoeff");
50
51 // Penalize negative parameters, just in case they show up
52 // when calculating the numeric derivative
53 if (H < std::numeric_limits<double>::epsilon() || D < std::numeric_limits<double>::epsilon()) {
54 for (size_t j = 0; j < nData; j++) {
55 out[j] = std::numeric_limits<double>::infinity();
56 }
57 return;
58 }
59
60 // Lorentzian intensities and HWHM
61 D *= 0.10; // conversion from 10^{-5}cm^2/s to Angstrom^2/ps, the internal
62 // units used
63 auto G = D * Q * Q;
64 for (size_t j = 0; j < nData; j++) {
65 auto E = xValues[j] - C;
66 out[j] += H * G / (G * G + E * E) / M_PI;
67 }
68}
69
76void FickDiffusionSQE::functionDeriv1D(Mantid::API::Jacobian *jacobian, const double *xValues, const size_t nData) {
77 const double deltaF{0.1}; // increase parameter by this fraction
78 const size_t nParam = this->nParams();
79 // cutoff defines the smallest change in the parameter when calculating the
80 // numerical derivative
81 std::map<std::string, double> cutoff;
82 cutoff["DiffCoeff"] = 0.2; // 0.2x10^(-5)cm^2/s
83 cutoff["Centre"] = 0.0001; // 0.1micro-eV
84 std::vector<double> out(nData);
85 this->applyTies();
86 this->function1D(out.data(), xValues, nData);
87
88 for (size_t iP = 0; iP < nParam; iP++) {
89 std::vector<double> derivative(nData);
90 if (this->isActive(iP)) {
91 const double pVal = this->activeParameter(iP);
92 const std::string pName = this->parameterName(iP);
93 if (pName == "Height") {
94 // exact derivative
95 this->setActiveParameter(iP, 1.0);
96 this->applyTies();
97 this->function1D(derivative.data(), xValues, nData);
98 } else {
99 // numerical derivative
100 double delta = cutoff[pName] > fabs(pVal * deltaF) ? cutoff[pName] : pVal * deltaF;
101 this->setActiveParameter(iP, pVal + delta);
102 this->applyTies();
103 this->function1D(derivative.data(), xValues, nData);
104 for (size_t i = 0; i < nData; i++) {
105 derivative[i] = (derivative[i] - out[i]) / delta;
106 }
107 }
108 this->setActiveParameter(iP, pVal); // restore the value of the parameter
109 // fill the jacobian for this parameter
110 for (size_t i = 0; i < nData; i++) {
111 jacobian->set(i, iP, derivative[i]);
112 }
113 }
114 }
115}
116
117} // namespace Mantid::CurveFitting::Functions
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
#define fabs(x)
Definition: Matrix.cpp:22
virtual void derivative(const FunctionDomain &domain, FunctionValues &values, const size_t order=1) const
Definition: IFunction1D.cpp:65
double asDouble() const
Returns double value if attribute is a double, throws exception otherwise.
Definition: IFunction.cpp:739
bool isActive(size_t i) const
Check if an active parameter i is actually active.
Definition: IFunction.cpp:160
virtual Attribute getAttribute(const std::string &name) const
Return a value of attribute attName.
Definition: IFunction.cpp:1394
virtual void applyTies()
Apply the ties.
Definition: IFunction.cpp:306
virtual double activeParameter(size_t i) const
Value of i-th active parameter.
Definition: IFunction.cpp:988
virtual void setActiveParameter(size_t i, double value)
Set new value of i-th active parameter.
Definition: IFunction.cpp:997
Represents the Jacobian in IFitFunction::functionDeriv.
Definition: Jacobian.h:22
virtual void set(size_t iY, size_t iP, double value)=0
Set a value to a Jacobian matrix element.
std::string parameterName(size_t i) const override
Returns the name of parameter i.
size_t nParams() const override
Total number of parameters.
Definition: ParamFunction.h:53
double getParameter(size_t i) const override
Get i-th parameter.
void function1D(double *out, const double *xValues, const size_t nData) const override
Calculate function values on an energy domain.
void functionDeriv1D(Mantid::API::Jacobian *jacobian, const double *xValues, const size_t nData) override
analytical/numerical derivative with respect to fitting parameters We carry out analytical derivative...
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
Kernel::Logger g_log("ExperimentInfo")
static logger object