Mantid
Loading...
Searching...
No Matches
Convolution.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//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
14#include "MantidAPI/IFunction.h"
17
18#include <algorithm>
19#include <cmath>
20#include <functional>
21
22#include <gsl/gsl_errno.h>
23#include <gsl/gsl_fft_halfcomplex.h>
24#include <gsl/gsl_fft_real.h>
25
26#include <fstream>
27#include <sstream>
28
29namespace {
30const double tolerance{0.02};
31}
32
34
35using namespace CurveFitting;
36using namespace Kernel;
37using namespace API;
38using std::placeholders::_1;
39
41 bool is1D = true;
42 const auto compositeFunction = std::dynamic_pointer_cast<CompositeFunction>(function);
43 if (compositeFunction) {
44 for (size_t index = 0; index < compositeFunction->nFunctions(); index++) {
45 if (!is1DCompositeFunction(compositeFunction->getFunction(index))) {
46 is1D = false;
47 break;
48 }
49 }
50 } else {
51 const IFunction1D_sptr fun = std::dynamic_pointer_cast<IFunction1D>(function);
52 is1D = fun ? true : false;
53 }
54 return is1D;
55}
56
57void evaluateFunctionOnRange(const IFunction_sptr &function, size_t domainSize, const double *range,
58 std::vector<double> &output) {
59 const IFunction1D_sptr fun = std::dynamic_pointer_cast<IFunction1D>(function);
60 if (!fun) {
61 FunctionDomain1DView ds(range, domainSize);
62 FunctionValues outs(ds);
63 function->function(ds, outs);
64 output = outs.toVector();
65 } else {
66 fun->function1D(output.data(), range, domainSize);
67 }
68}
69
71
72
74 declareAttribute("FixResolution", Attribute(true));
75 setAttributeValue("NumDeriv", true);
76}
77
79
80void Convolution::functionDeriv(const FunctionDomain &domain, Jacobian &jacobian) {
81 calNumericalDeriv(domain, jacobian);
82}
83
84void Convolution::setAttribute(const std::string &attName, const IFunction::Attribute &att) {
85 // if the resolution is there fix/unfix its parameters according to the
86 // attribute value
87 if (attName == "FixResolution" && nFunctions() > 0) {
88 bool fixRes = att.asBool();
89 auto &f = *getFunction(0);
90 for (size_t i = 0; i < f.nParams(); i++) {
91 if (fixRes) {
92 f.fix(i);
93 } else {
94 f.unfix(i);
95 }
96 }
97 }
99}
100
101namespace {
102// anonymous namespace for local definitions
103
104// A struct incapsulating workspaces for real fft
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);
111 }
112 gsl_fft_real_workspace *workspace;
113 gsl_fft_real_wavetable *wavetable;
114};
115} // namespace
116
125void Convolution::function(const FunctionDomain &domain, FunctionValues &values) const {
127
128 if (nFunctions() == 0) {
129 values.zeroCalculated();
130 return;
131 }
132 const auto &d1d = dynamic_cast<const FunctionDomain1D &>(domain);
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));
136 // positive x-values:
137 auto ixP = static_cast<size_t>(xValues[nData - 1] / dx);
138 auto ixN = nData - ixP - 1; // negative x-values (ixP+ixN=nData-1)
139
140 // determine whether to use FFT or Direct calculations
141 int assymmetry = abs(static_cast<int>(ixP - ixN));
142 if (xValues[0] * xValues[nData - 1] < 0 && assymmetry > tolerance * static_cast<double>(ixP + ixN)) {
143 functionDirectMode(domain, values);
144 } else {
145 functionFFTMode(domain, values);
146 }
147}
148
153 for (size_t index = 0; index < nFunctions(); index++) {
155 throw std::runtime_error("Convolution only enabled for 1D functions.");
156 }
157 }
158}
159
168 const auto &d1d = dynamic_cast<const FunctionDomain1D &>(domain);
169 size_t nData = domain.size();
170 const double *xValues = d1d.getPointerAt(0);
172 RealFFTWorkspace workspace(nData);
173 int n2 = static_cast<int>(nData) / 2;
174 bool odd = n2 * 2 != static_cast<int>(nData);
175 if (m_resolution.empty()) {
176 m_resolution.resize(nData);
177 // the resolution must be defined on interval -L < xr < L, L ==
178 // (xValues[nData-1] - xValues[0]) / 2
179 std::vector<double> xr(nData);
180 double dx = (xValues[nData - 1] - xValues[0]) / static_cast<double>((nData - 1));
181 // make sure that xr[nData/2] == 0.0
182 xr[n2] = 0.0;
183 for (int i = 1; i < n2; i++) {
184 double x = i * dx;
185 xr[n2 + i] = x;
186 xr[n2 - i] = -x;
187 }
188
189 xr[0] = -n2 * dx;
190 if (odd)
191 xr[nData - 1] = -xr[0];
193
194 // rotate the data to produce the right transform
195 if (odd) {
196 double tmp = m_resolution[nData - 1];
197 for (int i = n2 - 1; i >= 0; i--) {
198 m_resolution[n2 + i + 1] = m_resolution[i];
199 m_resolution[i] = m_resolution[n2 + i];
200 }
201 m_resolution[n2] = tmp;
202 } else {
203 for (int i = 0; i < n2; i++) {
204 double tmp = m_resolution[i];
205 m_resolution[i] = m_resolution[n2 + i];
206 m_resolution[n2 + i] = tmp;
207 }
208 }
209 gsl_fft_real_transform(m_resolution.data(), 1, nData, workspace.wavetable, workspace.workspace);
210 std::transform(m_resolution.begin(), m_resolution.end(), m_resolution.begin(),
211 std::bind(std::multiplies<double>(), _1, dx));
212 }
213
214 // Now m_resolution contains fourier transform of the resolution
215
216 if (nFunctions() == 1) {
217 // return the resolution transform for testing
218 double dx = 1.; // nData > 1? xValues[1] - xValues[0]: 1.;
219 std::transform(m_resolution.begin(), m_resolution.end(), values.getPointerToCalculated(0),
220 std::bind(std::multiplies<double>(), _1, dx));
221 return;
222 }
223
224 // check for delta functions
225 std::vector<std::shared_ptr<DeltaFunction>> dltFuns;
226 double dltF = 0;
227 bool deltaFunctionsOnly = false;
228 bool deltaShifted = false;
229 CompositeFunction_sptr cf = std::dynamic_pointer_cast<CompositeFunction>(getFunction(1));
230 if (cf) {
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));
234 if (df) {
235 dltFuns.emplace_back(df);
236 if (df->getParameter("Centre") != 0.0) {
237 deltaShifted = true;
238 }
239 dltF += df->getParameter("Height") * df->HeightPrefactor();
240 }
241 }
242 if (dltFuns.size() == cf->nFunctions()) {
243 // all delta functions - return scaled resolution
244 deltaFunctionsOnly = true;
245 }
246 } else if (auto df = std::dynamic_pointer_cast<DeltaFunction>(getFunction(1))) {
247 // single delta function - return scaled resolution
248 deltaFunctionsOnly = true;
249 dltFuns.emplace_back(df);
250 if (df->getParameter("Centre") != 0.0) {
251 deltaShifted = true;
252 }
253 dltF = df->getParameter("Height") * df->HeightPrefactor();
254 }
255
256 // out points to the calculated values in values
257 double *out = values.getPointerToCalculated(0);
258
259 if (!deltaFunctionsOnly) {
260 // Transform the model function
261 getFunction(1)->function(domain, values);
262 gsl_fft_real_transform(out, 1, nData, workspace.wavetable, workspace.workspace);
263
264 // Fourier transform is integration - multiply by the step in the
265 // integration variable
266 double dx = nData > 1 ? xValues[1] - xValues[0] : 1.;
267 std::transform(out, out + nData, out, std::bind(std::multiplies<double>(), _1, dx));
268
269 // now out contains fourier transform of the model function
270
271 HalfComplex res(m_resolution.data(), nData);
272 HalfComplex fun(out, nData);
273
274 // Multiply transforms of the resolution and model functions
275 // Result is stored in fun
276 for (size_t i = 0; i <= res.size(); i++) {
277 // complex multiplication
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);
283 }
284
285 // Inverse fourier transform of fun
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);
289
290 // Inverse fourier transform is integration - multiply by the step in the
291 // integration variable
292 dx = nData > 1 ? 1. / (xValues[1] - xValues[0]) : 1.;
293 std::transform(out, out + nData, out, std::bind(std::multiplies<double>(), _1, dx));
294 } else {
295 values.zeroCalculated();
296 }
297
298 if (dltF != 0.0 && !deltaShifted) {
299 // If model contains any delta functions their effect is addition of scaled
300 // resolution
301 std::vector<double> tmp(nData);
302 evaluateFunctionOnRange(getFunction(0), nData, xValues, tmp);
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);
312 evaluateFunctionOnRange(getFunction(0), nData, &x[0], tmp);
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>());
315 }
316 }
317
318} // end of Convolution::functionFFTMode
319
330 const auto &d1d = dynamic_cast<const FunctionDomain1D &>(domain);
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); // positive
335 // x-values
336 auto ixN = nData - ixP - 1; // negative x-values (ixP+ixN=nData-1)
337
339
340 // double the domain where to evaluate the convolution. Guarantees complete
341 // overlap betwen convolution and signal in the original range.
342 const size_t mData = nData + ixN + ixP; // equal to 2*nData-1
343 std::vector<double> xValuesExtdVec(mData);
344 double *xValuesExtd = xValuesExtdVec.data();
345 Mantid::API::FunctionDomain1DView domainExtd(xValuesExtd, mData);
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;
349 }
350
351 if (m_resolution.empty()) {
352 m_resolution.resize(nData);
353 }
354 // Fill m_resolution with the resolution function data
355 // Lines 341-349 is duplicated in functionFFTmode. To be cleanup
356 // in issue 16064
357 evaluateFunctionOnRange(getFunction(0), nData, &xValues[0], m_resolution);
358
359 // Reverse the axis of the resolution data
360 std::reverse(m_resolution.begin(), m_resolution.end());
361
362 // check for delta functions
363 std::vector<std::shared_ptr<DeltaFunction>> dltFuns;
364 double dltF = 0;
365 bool deltaFunctionsOnly = false;
366 bool deltaShifted = false;
367 CompositeFunction_sptr cf = std::dynamic_pointer_cast<CompositeFunction>(getFunction(1));
368 if (cf) {
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));
372 if (df) {
373 dltFuns.emplace_back(df);
374 if (df->getParameter("Centre") != 0.0) {
375 deltaShifted = true;
376 }
377 dltF += df->getParameter("Height") * df->HeightPrefactor();
378 }
379 }
380 if (dltFuns.size() == cf->nFunctions()) {
381 // all delta functions - return scaled resolution
382 deltaFunctionsOnly = true;
383 }
384 } else if (auto df = std::dynamic_pointer_cast<DeltaFunction>(getFunction(1))) {
385 // single delta function - return scaled resolution
386 deltaFunctionsOnly = true;
387 dltFuns.emplace_back(df);
388 if (df->getParameter("Centre") != 0.0) {
389 deltaShifted = true;
390 }
391 dltF = df->getParameter("Height") * df->HeightPrefactor();
392 }
393
394 // store the result of the convolution
395 double *out = values.getPointerToCalculated(0);
396
397 if (!deltaFunctionsOnly) {
398 // Evaluate the model on the extended domain
399 Mantid::API::FunctionValues valuesExtd(domainExtd);
400 getFunction(1)->function(domainExtd, valuesExtd);
401 // Convolve with resolution
402 double *outExt = valuesExtd.getPointerToCalculated(0);
403 for (size_t i = 0; i < nData; i++) {
404 double tmp{0.0};
405 for (size_t j = 0; j < nData; j++) {
406 tmp += outExt[i + j] * m_resolution[j];
407 }
408 out[i] = tmp * dx;
409 }
410 } else {
411 values.zeroCalculated();
412 }
413
414 if (dltF != 0.0 && !deltaShifted) {
415 // If model contains any delta functions their effect is addition of scaled
416 // resolution
417 // Lines 412-430 is duplicated in functionFFTmode. To be cleanup
418 // in issue 16064
419 std::vector<double> tmp(nData);
420 evaluateFunctionOnRange(getFunction(0), nData, xValues, tmp);
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);
430 evaluateFunctionOnRange(getFunction(0), nData, &x[0], tmp);
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>());
433 }
434 }
435
436} // end of Convolution::functionDirectMode()
437
448 if (std::dynamic_pointer_cast<Convolution>(f)) {
449 throw std::runtime_error("Nested convolutions are not allowed.");
450 }
451 if (nFunctions() == 0 && getAttribute("FixResolution").asBool()) {
452 for (size_t i = 0; i < f->nParams(); i++) {
453 f->fix(i);
454 }
455 }
456 size_t iFun = 0;
457 if (nFunctions() < 2) {
459 } else {
461 if (!f1) {
462 throw std::runtime_error("IFunction expected but function of another type found");
463 }
464 CompositeFunction_sptr cf = std::dynamic_pointer_cast<CompositeFunction>(f1);
465 if (cf == nullptr) {
466 cf = std::dynamic_pointer_cast<CompositeFunction>(
467 API::FunctionFactory::Instance().createFunction("CompositeFunction"));
469 cf->addFunction(f1);
471 }
472 cf->addFunction(f);
474 iFun = 1;
475 }
476 return iFun;
477}
478
484
488 // refresh when calculation for the first time
489 bool needRefreshing = m_resolution.empty();
490 if (!needRefreshing) {
491 // if resolution has active parameters always refresh
492 IFunction &res = *getFunction(0);
493 for (size_t i = 0; i < res.nParams(); ++i) {
494 if (res.isActive(i)) {
495 needRefreshing = true;
496 break;
497 }
498 }
499 }
500 if (!needRefreshing)
501 return;
502 // delete fourier transform of the resolution to force its recalculation
503 m_resolution.clear();
504}
505
506} // namespace Mantid::CurveFitting::Functions
gsl_vector * tmp
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
Definition: IndexPeaks.cpp:114
Mantid::API::IFunction::Attribute Attribute
Definition: IsoRotDiff.cpp:25
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
double tolerance
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.
Definition: IFunction.h:282
bool asBool() const
Returns bool value if attribute is a bool, throws exception otherwise.
Definition: IFunction.cpp:752
This is an interface to a fitting function - a semi-abstarct class.
Definition: IFunction.h:163
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.
Definition: IFunction.cpp:160
void calNumericalDeriv(const FunctionDomain &domain, Jacobian &jacobian)
Calculate numerical derivatives.
Definition: IFunction.cpp:1031
Represents the Jacobian in IFitFunction::functionDeriv.
Definition: Jacobian.h:22
Class for helping to read the transformed data.
Definition: Convolution.h:39
size_t size() const
Returns the size of the transform.
Definition: Convolution.h:51
double real(size_t i) const
The real part of i-th transform coefficient.
Definition: Convolution.h:57
void set(size_t i, const double &re, const double &im)
Set a new value for i-th complex coefficient.
Definition: Convolution.h:84
double imag(size_t i) const
The imaginary part of i-th transform coefficient.
Definition: Convolution.h:69
Performes convolution of two functions.
Definition: Convolution.h:28
void functionDeriv(const API::FunctionDomain &domain, API::Jacobian &jacobian) override
Derivatives of function with respect to active parameters.
Definition: Convolution.cpp:80
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
Definition: Convolution.cpp:78
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.
Definition: Convolution.cpp:84
std::vector< double > m_resolution
Keep the Fourier transform of the resolution function (divided by the step in xValues) when in FFT mo...
Definition: Convolution.h:135
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
Definition: IFunction1D.h:80
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
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)
Definition: Convolution.cpp:57
bool is1DCompositeFunction(const IFunction_sptr &function)
Definition: Convolution.cpp:40