Mantid
Loading...
Searching...
No Matches
FitMW.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// Includes
8//----------------------------------------------------------------------
14
22#include "MantidAPI/TextAxis.h"
24
29#include "MantidKernel/Matrix.h"
30
31#include <algorithm>
32#include <cmath>
33#include <numeric>
34
35namespace Mantid::CurveFitting {
36
37using namespace Kernel;
38using API::MatrixWorkspace;
39using API::Workspace;
40
41namespace {
42
46struct RangePoint {
47 enum Kind : char { Openning, Closing };
49 Kind kind;
51 double value;
54 bool operator<(const RangePoint &point) const {
55 if (this->value == point.value) {
56 // If an Openning and Closing points have the same value
57 // the Openning one should go first (be "smaller").
58 // This way the procedure of joinOverlappingRanges will join
59 // the ranges that meet at these points.
60 return this->kind == Openning;
61 }
62 return this->value < point.value;
63 }
64};
65
70void joinOverlappingRanges(std::vector<double> &exclude) {
71 if (exclude.empty()) {
72 return;
73 }
74 // The situation here is similar to matching brackets in an expression.
75 // If we sort all the points in the input vector remembering their kind
76 // then a separate exclusion region starts with the first openning bracket
77 // and ends with the matching closing bracket. All brackets (points) inside
78 // them can be dropped.
79
80 // Wrap the points into helper struct RangePoint
81 std::vector<RangePoint> points;
82 points.reserve(exclude.size());
83 for (auto point = exclude.begin(); point != exclude.end(); point += 2) {
84 points.emplace_back(RangePoint{RangePoint::Openning, *point});
85 points.emplace_back(RangePoint{RangePoint::Closing, *(point + 1)});
86 }
87 // Sort the points according to the operator defined in RangePoint.
88 std::sort(points.begin(), points.end());
89
90 // Clear the argument vector.
91 exclude.clear();
92 // Start the level counter which shows the number of unmatched openning
93 // brackets.
94 size_t level(0);
95 for (auto &point : points) {
96 if (point.kind == RangePoint::Openning) {
97 if (level == 0) {
98 // First openning bracket starts a new exclusion range.
99 exclude.emplace_back(point.value);
100 }
101 // Each openning bracket increases the level
102 ++level;
103 } else {
104 if (level == 1) {
105 // The bracket that makes level 0 is an end of a range.
106 exclude.emplace_back(point.value);
107 }
108 // Each closing bracket decreases the level
109 --level;
110 }
111 }
112}
113
114} // namespace
115
123FitMW::FitMW(Kernel::IPropertyManager *fit, const std::string &workspacePropertyName, FitMW::DomainType domainType)
124 : IMWDomainCreator(fit, workspacePropertyName, domainType), m_maxSize(0), m_normalise(false) {}
125
133 : IMWDomainCreator(nullptr, std::string(), domainType), m_maxSize(10), m_normalise(false) {}
134
142 // if property manager is set overwrite any set parameters
143 if (m_manager) {
144
145 if (m_domainType != Simple) {
146 const int maxSizeInt = m_manager->getProperty(m_maxSizePropertyName);
147 m_maxSize = static_cast<size_t>(maxSizeInt);
148 }
151 if (m_exclude.size() % 2 != 0) {
152 throw std::runtime_error("Exclude property has an odd number of entries. "
153 "It has to be even as each pair specifies a "
154 "start and an end of an interval to exclude.");
155 }
156 joinOverlappingRanges(m_exclude);
157 }
158}
159
166void FitMW::declareDatasetProperties(const std::string &suffix, bool addProp) {
168 m_maxSizePropertyName = "MaxSize" + suffix;
169 m_normalisePropertyName = "Normalise" + suffix;
170 m_excludePropertyName = "Exclude" + suffix;
171
172 if (addProp) {
174 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
175 mustBePositive->setLower(0);
177 "The maximum number of values per a simple domain.");
178 }
181 "An option to normalise the histogram data (divide be the bin "
182 "width).");
183 }
185 auto mustBeOrderedPairs = std::make_shared<ArrayOrderedPairsValidator<double>>();
187 "A list of pairs of doubles that specify ranges that "
188 "must be excluded from fit.");
189 }
190 }
191}
192
194void FitMW::createDomain(std::shared_ptr<API::FunctionDomain> &domain, std::shared_ptr<API::FunctionValues> &values,
195 size_t i0) {
197
198 const auto &X = m_matrixWorkspace->x(m_workspaceIndex);
199
200 // find the fitting interval: from -> to
201 size_t endIndex = 0;
202 std::tie(m_startIndex, endIndex) = getXInterval();
203 auto from = X.begin() + m_startIndex;
204 auto to = X.begin() + endIndex;
205 auto n = endIndex - m_startIndex;
206
207 if (m_domainType != Simple) {
208 if (m_maxSize < n) {
210 domain.reset(seqDomain);
211 size_t m = 0;
212 while (m < n) {
213 // create a simple creator
214 auto creator = new FitMW;
216 creator->setWorkspaceIndex(m_workspaceIndex);
217 size_t k = m + m_maxSize;
218 if (k > n)
219 k = n;
220 creator->setRange(*(from + m), *(from + k - 1));
221 seqDomain->addCreator(API::IDomainCreator_sptr(creator));
222 m = k;
223 }
224 values.reset();
225 return;
226 }
227 // else continue with simple domain
228 }
229
230 // set function domain
231 if (m_matrixWorkspace->isHistogramData()) {
232 std::vector<double> x(static_cast<size_t>(to - from));
233 auto it = from;
234 for (size_t i = 0; it != to; ++it, ++i) {
235 x[i] = (*it + *(it + 1)) / 2;
236 }
238 x.clear();
239 } else {
240 domain.reset(new API::FunctionDomain1DSpectrum(m_workspaceIndex, from, to));
241 }
242
243 if (!values) {
244 values.reset(new API::FunctionValues(*domain));
245 } else {
246 values->expand(i0 + domain->size());
247 }
248
249 bool shouldNormalise = m_normalise && m_matrixWorkspace->isHistogramData();
250
251 // set the data to fit to
252 assert(n == domain->size());
253 const auto &Y = m_matrixWorkspace->y(m_workspaceIndex);
254 const auto &E = m_matrixWorkspace->e(m_workspaceIndex);
255 if (endIndex > Y.size()) {
256 throw std::runtime_error("FitMW: Inconsistent MatrixWorkspace");
257 }
258
259 // Helps find points excluded form fit.
260 ExcludeRangeFinder excludeFinder(m_exclude, X.front(), X.back());
261
262 for (size_t i = m_startIndex; i < endIndex; ++i) {
263 size_t j = i - m_startIndex + i0;
264 double y = Y[i];
265 double error = E[i];
266 double weight = 0.0;
267
268 if (shouldNormalise) {
269 double binWidth = X[i + 1] - X[i];
270 if (binWidth == 0.0) {
271 throw std::runtime_error("Zero width bin found, division by zero.");
272 }
273 y /= binWidth;
274 error /= binWidth;
275 }
276
277 if (excludeFinder.isExcluded(X[i])) {
278 weight = 0.0;
279 } else if (!std::isfinite(y)) {
280 // nan or inf data
282 throw std::runtime_error("Infinte number or NaN found in input data.");
283 y = 0.0; // leaving inf or nan would break the fit
284 } else if (!std::isfinite(error)) {
285 // nan or inf error
287 throw std::runtime_error("Infinte number or NaN found in input data.");
288 } else if (error <= 0) {
290 weight = 1.0;
291 } else {
292 weight = 1.0 / error;
293 if (!std::isfinite(weight)) {
295 throw std::runtime_error("Error of a data point is probably too small.");
296 weight = 0.0;
297 }
298 }
299
300 values->setFitData(j, y);
301 values->setFitWeight(j, weight);
302 }
303 m_domain = std::dynamic_pointer_cast<API::FunctionDomain1D>(domain);
304 m_values = values;
305}
306
315std::shared_ptr<API::Workspace> FitMW::createOutputWorkspace(const std::string &baseName, API::IFunction_sptr function,
316 std::shared_ptr<API::FunctionDomain> domain,
317 std::shared_ptr<API::FunctionValues> values,
318 const std::string &outputWorkspacePropertyName) {
319 auto ws = IMWDomainCreator::createOutputWorkspace(baseName, function, domain, values, outputWorkspacePropertyName);
320 auto &mws = dynamic_cast<MatrixWorkspace &>(*ws);
321
322 if (m_normalise && m_matrixWorkspace->isHistogramData()) {
323 const auto &X = mws.x(0);
324 std::vector<double> binWidths(X.size());
325 std::adjacent_difference(X.begin(), X.end(), binWidths.begin());
326 for (size_t ispec = 1; ispec < mws.getNumberHistograms(); ++ispec) {
327 auto &Y = mws.mutableY(ispec);
328 std::transform(binWidths.begin() + 1, binWidths.end(), Y.begin(), Y.begin(), std::multiplies<double>());
329 auto &E = mws.mutableE(ispec);
330 std::transform(binWidths.begin() + 1, binWidths.end(), E.begin(), E.begin(), std::multiplies<double>());
331 }
332 }
333 return ws;
334}
335
336} // namespace Mantid::CurveFitting
Kind kind
The kind of the point: either openning or closing the range.
Definition: FitMW.cpp:49
double value
The value of the point.
Definition: FitMW.cpp:51
double error
Definition: IndexPeaks.cpp:133
Specialization of FunctionDomain1DVector for spectra of MatrixWorkspaces.
A class to store values calculated by a function.
DomainType
Type of domain to create.
Kernel::IPropertyManager * m_manager
Pointer to a property manager.
bool m_ignoreInvalidData
Flag to ignore nans, infinities and zero errors.
void declareProperty(Kernel::Property *prop, const std::string &doc)
Declare a property to the algorithm.
DomainType m_domainType
Domain type.
Base MatrixWorkspace Abstract Class.
ExcludeRangeFinder : Helper clss that finds if a point should be excluded from fit.
bool isExcluded(double value)
Check if an x-value lies in an exclusion range.
Creates FunctionDomain1D form a spectrum in a MatrixWorkspace.
Definition: FitMW.h:38
std::string m_excludePropertyName
Store the Exclude property name.
Definition: FitMW.h:74
bool m_normalise
Option to normalise the data.
Definition: FitMW.h:79
void createDomain(std::shared_ptr< API::FunctionDomain > &domain, std::shared_ptr< API::FunctionValues > &values, size_t i0=0) override
Create a domain from the input workspace.
Definition: FitMW.cpp:194
void setParameters() const override
Set all parameters.
Definition: FitMW.cpp:140
std::vector< double > m_exclude
Ranges that must be excluded from fit.
Definition: FitMW.h:81
std::string m_maxSizePropertyName
Store maxSize property name.
Definition: FitMW.h:70
size_t m_maxSize
Max size for seq domain.
Definition: FitMW.h:77
void declareDatasetProperties(const std::string &suffix="", bool addProp=true) override
Declare properties that specify the dataset within the workspace to fit to.
Definition: FitMW.cpp:166
FitMW(Kernel::IPropertyManager *fit, const std::string &workspacePropertyName, DomainType domainType=Simple)
Constructor.
Definition: FitMW.cpp:123
std::string m_normalisePropertyName
Store Normalise property name.
Definition: FitMW.h:72
std::shared_ptr< API::Workspace > createOutputWorkspace(const std::string &baseName, API::IFunction_sptr function, std::shared_ptr< API::FunctionDomain > domain, std::shared_ptr< API::FunctionValues > values, const std::string &outputWorkspacePropertyName) override
Create an output workspace.
Definition: FitMW.cpp:315
A base class for domain creators taking 1D data from a spectrum of a matrix workspace.
void declareDatasetProperties(const std::string &suffix="", bool addProp=true) override
Declare properties that specify the dataset within the workspace to fit to.
void setWorkspace(std::shared_ptr< API::MatrixWorkspace > ws)
Set the workspace.
std::weak_ptr< API::FunctionDomain1D > m_domain
Store the created domain and values.
std::pair< size_t, size_t > getXInterval() const
Calculate size and starting iterator in the X array.
std::shared_ptr< API::MatrixWorkspace > m_matrixWorkspace
The input MareixWorkspace.
std::weak_ptr< API::FunctionValues > m_values
virtual void setParameters() const
Set all parameters.
size_t m_workspaceIndex
The workspace index.
std::shared_ptr< API::Workspace > createOutputWorkspace(const std::string &baseName, API::IFunction_sptr function, std::shared_ptr< API::FunctionDomain > domain, std::shared_ptr< API::FunctionValues > values, const std::string &outputWorkspacePropertyName) override
Create an output workspace.
An implementation of CompositeDomain.
Definition: SeqDomain.h:30
void addCreator(const API::IDomainCreator_sptr &creator)
Add new domain creator.
Definition: SeqDomain.cpp:48
static SeqDomain * create(API::IDomainCreator::DomainType type)
Create an instance of SeqDomain in one of two forms: either SeqDomain for sequential domain creation ...
Definition: SeqDomain.cpp:60
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
Interface to PropertyManager.
virtual bool existsProperty(const std::string &name) const =0
Checks whether the named property is already in the list of managed property.
virtual TypedValue getProperty(const std::string &name) const =0
Get the value of a property.
The concrete, templated class for properties.
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
std::shared_ptr< IDomainCreator > IDomainCreator_sptr
Typedef for a shared pointer to IDomainCreator.
STL namespace.
constexpr bool operator<(const wide_integer< Bits, Signed > &lhs, const wide_integer< Bits2, Signed2 > &rhs)