22#include <gsl/gsl_errno.h>
23#include <gsl/gsl_fft_halfcomplex.h>
24#include <gsl/gsl_fft_real.h>
35using namespace CurveFitting;
36using namespace Kernel;
38using std::placeholders::_1;
42 const auto compositeFunction = std::dynamic_pointer_cast<CompositeFunction>(function);
43 if (compositeFunction) {
44 for (
size_t index = 0;
index < compositeFunction->nFunctions();
index++) {
51 const IFunction1D_sptr fun = std::dynamic_pointer_cast<IFunction1D>(function);
52 is1D = fun ? true :
false;
58 std::vector<double> &output) {
59 const IFunction1D_sptr fun = std::dynamic_pointer_cast<IFunction1D>(function);
63 function->function(ds, outs);
66 fun->function1D(output.data(), range, domainSize);
74 declareAttribute(
"FixResolution",
Attribute(
true));
75 setAttributeValue(
"NumDeriv",
true);
87 if (attName ==
"FixResolution" &&
nFunctions() > 0) {
88 bool fixRes = att.
asBool();
90 for (
size_t i = 0; i < f.nParams(); i++) {
105struct RealFFTWorkspace {
106 explicit RealFFTWorkspace(
size_t nData)
107 :
workspace(gsl_fft_real_workspace_alloc(nData)),
wavetable(gsl_fft_real_wavetable_alloc(nData)) {}
108 ~RealFFTWorkspace() {
109 gsl_fft_real_wavetable_free(wavetable);
110 gsl_fft_real_workspace_free(workspace);
133 const size_t nData = domain.
size();
134 const double *xValues = d1d.getPointerAt(0);
135 double dx = (xValues[nData - 1] - xValues[0]) /
static_cast<double>((nData - 1));
137 auto ixP =
static_cast<size_t>(xValues[nData - 1] / dx);
138 auto ixN = nData - ixP - 1;
141 int assymmetry = abs(
static_cast<int>(ixP - ixN));
142 if (xValues[0] * xValues[nData - 1] < 0 && assymmetry >
tolerance *
static_cast<double>(ixP + ixN)) {
155 throw std::runtime_error(
"Convolution only enabled for 1D functions.");
169 size_t nData = domain.
size();
170 const double *xValues = d1d.getPointerAt(0);
173 int n2 =
static_cast<int>(nData) / 2;
174 bool odd = n2 * 2 !=
static_cast<int>(nData);
179 std::vector<double> xr(nData);
180 double dx = (xValues[nData - 1] - xValues[0]) /
static_cast<double>((nData - 1));
183 for (
int i = 1; i < n2; i++) {
191 xr[nData - 1] = -xr[0];
197 for (
int i = n2 - 1; i >= 0; i--) {
203 for (
int i = 0; i < n2; i++) {
211 std::bind(std::multiplies<double>(), _1, dx));
220 std::bind(std::multiplies<double>(), _1, dx));
225 std::vector<std::shared_ptr<DeltaFunction>> dltFuns;
227 bool deltaFunctionsOnly =
false;
228 bool deltaShifted =
false;
231 dltFuns.reserve(cf->nFunctions());
232 for (
size_t i = 0; i < cf->nFunctions(); ++i) {
233 auto df = std::dynamic_pointer_cast<DeltaFunction>(cf->getFunction(i));
235 dltFuns.emplace_back(df);
236 if (df->getParameter(
"Centre") != 0.0) {
239 dltF += df->getParameter(
"Height") * df->HeightPrefactor();
242 if (dltFuns.size() == cf->nFunctions()) {
244 deltaFunctionsOnly =
true;
246 }
else if (
auto df = std::dynamic_pointer_cast<DeltaFunction>(
getFunction(1))) {
248 deltaFunctionsOnly =
true;
249 dltFuns.emplace_back(df);
250 if (df->getParameter(
"Centre") != 0.0) {
253 dltF = df->getParameter(
"Height") * df->HeightPrefactor();
259 if (!deltaFunctionsOnly) {
266 double dx = nData > 1 ? xValues[1] - xValues[0] : 1.;
267 std::transform(out, out + nData, out, std::bind(std::multiplies<double>(), _1, dx));
276 for (
size_t i = 0; i <= res.
size(); i++) {
278 double res_r = res.
real(i);
279 double res_i = res.
imag(i);
280 double fun_r = fun.
real(i);
281 double fun_i = fun.
imag(i);
282 fun.
set(i, res_r * fun_r - res_i * fun_i, res_r * fun_i + res_i * fun_r);
286 gsl_fft_halfcomplex_wavetable *wavetable_r = gsl_fft_halfcomplex_wavetable_alloc(nData);
287 gsl_fft_halfcomplex_inverse(out, 1, nData, wavetable_r,
workspace.workspace);
288 gsl_fft_halfcomplex_wavetable_free(wavetable_r);
292 dx = nData > 1 ? 1. / (xValues[1] - xValues[0]) : 1.;
293 std::transform(out, out + nData, out, std::bind(std::multiplies<double>(), _1, dx));
298 if (dltF != 0.0 && !deltaShifted) {
301 std::vector<double>
tmp(nData);
303 std::transform(
tmp.begin(),
tmp.end(),
tmp.begin(), std::bind(std::multiplies<double>(), _1, dltF));
304 std::transform(out, out + nData,
tmp.data(), out, std::plus<double>());
305 }
else if (!dltFuns.empty()) {
306 std::vector<double>
x(nData);
307 for (
const auto &df : dltFuns) {
308 double shift = -df->getParameter(
"Centre");
309 dltF = df->getParameter(
"Height") * df->HeightPrefactor();
310 std::transform(xValues, xValues + nData,
x.data(), std::bind(std::plus<double>(), _1, shift));
311 std::vector<double>
tmp(nData);
313 std::transform(
tmp.begin(),
tmp.end(),
tmp.begin(), std::bind(std::multiplies<double>(), _1, dltF));
314 std::transform(out, out + nData,
tmp.data(), out, std::plus<double>());
331 const size_t nData = domain.
size();
332 const double *xValues = d1d.getPointerAt(0);
333 double dx = (xValues[nData - 1] - xValues[0]) /
static_cast<double>((nData - 1));
334 auto ixP =
static_cast<size_t>(xValues[nData - 1] / dx);
336 auto ixN = nData - ixP - 1;
342 const size_t mData = nData + ixN + ixP;
343 std::vector<double> xValuesExtdVec(mData);
344 double *xValuesExtd = xValuesExtdVec.data();
346 double Dx = dx *
static_cast<double>(ixN + ixP);
347 for (
size_t i = 0; i < mData; i++) {
348 xValuesExtd[i] = -Dx +
static_cast<double>(i) * dx;
363 std::vector<std::shared_ptr<DeltaFunction>> dltFuns;
365 bool deltaFunctionsOnly =
false;
366 bool deltaShifted =
false;
369 dltFuns.reserve(cf->nFunctions());
370 for (
size_t i = 0; i < cf->nFunctions(); ++i) {
371 auto df = std::dynamic_pointer_cast<DeltaFunction>(cf->getFunction(i));
373 dltFuns.emplace_back(df);
374 if (df->getParameter(
"Centre") != 0.0) {
377 dltF += df->getParameter(
"Height") * df->HeightPrefactor();
380 if (dltFuns.size() == cf->nFunctions()) {
382 deltaFunctionsOnly =
true;
384 }
else if (
auto df = std::dynamic_pointer_cast<DeltaFunction>(
getFunction(1))) {
386 deltaFunctionsOnly =
true;
387 dltFuns.emplace_back(df);
388 if (df->getParameter(
"Centre") != 0.0) {
391 dltF = df->getParameter(
"Height") * df->HeightPrefactor();
397 if (!deltaFunctionsOnly) {
403 for (
size_t i = 0; i < nData; i++) {
405 for (
size_t j = 0; j < nData; j++) {
414 if (dltF != 0.0 && !deltaShifted) {
419 std::vector<double>
tmp(nData);
421 std::transform(
tmp.begin(),
tmp.end(),
tmp.begin(), std::bind(std::multiplies<double>(), _1, dltF));
422 std::transform(out, out + nData,
tmp.data(), out, std::plus<double>());
423 }
else if (!dltFuns.empty()) {
424 std::vector<double>
x(nData);
425 for (
const auto &df : dltFuns) {
426 double shift = -df->getParameter(
"Centre");
427 dltF = df->getParameter(
"Height") * df->HeightPrefactor();
428 std::transform(xValues, xValues + nData,
x.data(), std::bind(std::plus<double>(), _1, shift));
429 std::vector<double>
tmp(nData);
431 std::transform(
tmp.begin(),
tmp.end(),
tmp.begin(), std::bind(std::multiplies<double>(), _1, dltF));
432 std::transform(out, out + nData,
tmp.data(), out, std::plus<double>());
448 if (std::dynamic_pointer_cast<Convolution>(f)) {
449 throw std::runtime_error(
"Nested convolutions are not allowed.");
452 for (
size_t i = 0; i < f->nParams(); i++) {
462 throw std::runtime_error(
"IFunction expected but function of another type found");
466 cf = std::dynamic_pointer_cast<CompositeFunction>(
490 if (!needRefreshing) {
493 for (
size_t i = 0; i < res.
nParams(); ++i) {
495 needRefreshing =
true;
gsl_fft_real_wavetable * wavetable
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
IPeaksWorkspace_sptr workspace
Mantid::API::IFunction::Attribute Attribute
std::map< DeltaEMode::Type, std::string > index
virtual size_t addFunction(IFunction_sptr f)
Add a function at the back of the internal function list.
std::size_t nFunctions() const override
Number of functions.
Attribute getAttribute(const std::string &name) const override
Return a value of attribute attName.
void removeFunction(size_t i)
Remove a function.
IFunction_sptr getFunction(std::size_t i) const override
Returns the pointer to i-th function.
void setAttribute(const std::string &name, const API::IFunction::Attribute &value) override
Set a value of a named attribute.
void checkFunction()
Check the function.
1D domain - a wrapper around an array of doubles.
Represent a domain for functions of one real argument.
Base class that represents the domain of a function.
virtual size_t size() const =0
Return the number of points in the domain.
A class to store values calculated by a function.
void zeroCalculated()
Set all calculated values to zero.
double * getPointerToCalculated(size_t i)
Get a pointer to calculated data at index i.
std::vector< double > toVector() const
Return the calculated values as a vector.
Attribute is a non-fitting parameter.
bool asBool() const
Returns bool value if attribute is a bool, throws exception otherwise.
This is an interface to a fitting function - a semi-abstarct class.
virtual size_t nParams() const =0
Total number of parameters.
bool isActive(size_t i) const
Check if an active parameter i is actually active.
void calNumericalDeriv(const FunctionDomain &domain, Jacobian &jacobian)
Calculate numerical derivatives.
Represents the Jacobian in IFitFunction::functionDeriv.
Class for helping to read the transformed data.
size_t size() const
Returns the size of the transform.
double real(size_t i) const
The real part of i-th transform coefficient.
void set(size_t i, const double &re, const double &im)
Set a new value for i-th complex coefficient.
double imag(size_t i) const
The imaginary part of i-th transform coefficient.
Performes convolution of two functions.
void functionDeriv(const API::FunctionDomain &domain, API::Jacobian &jacobian) override
Derivatives of function with respect to active parameters.
void functionFFTMode(const API::FunctionDomain &domain, API::FunctionValues &values) const
Calculates convolution of the two member functions when the domain is symmetric with respect to inver...
size_t addFunction(API::IFunction_sptr f) override
Add a function.
void innerFunctionsAre1D() const
Checks that the functions being convolted are 1D and throws if not.
void init() override
overwrite IFunction base class method, which declare function parameters
void refreshResolution() const
Deletes and zeroes pointer m_resolution forsing function(...) to recalculate the resolution function.
void function(const API::FunctionDomain &domain, API::FunctionValues &values) const override
Function you want to fit to.
void setUpForFit() override
Set up the function for a fit.
void functionDirectMode(const API::FunctionDomain &domain, API::FunctionValues &values) const
Calculates convolution of the two member functions when the domain is not symmetric with respect to i...
void setAttribute(const std::string &attName, const Attribute &) override
Set a value to attribute attName.
std::vector< double > m_resolution
Keep the Fourier transform of the resolution function (divided by the step in xValues) when in FFT mo...
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< IFunction1D > IFunction1D_sptr
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
std::shared_ptr< CompositeFunction > CompositeFunction_sptr
shared pointer to the composite function base class
double f1(const double x, const double G, const double w0)
void evaluateFunctionOnRange(const IFunction_sptr &function, size_t domainSize, const double *range, std::vector< double > &output)
bool is1DCompositeFunction(const IFunction_sptr &function)