Mantid
Loading...
Searching...
No Matches
CalculatePolynomialBackground.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 +
8
19
20#include <utility>
21
22namespace {
24namespace Prop {
25const static std::string COST_FUNCTION = "CostFunction";
26const static std::string INPUT_WS = "InputWorkspace";
27const static std::string OUTPUT_WS = "OutputWorkspace";
28const static std::string POLY_DEGREE = "Degree";
29const static std::string XRANGES = "XRanges";
30const static std::string MINIMIZER = "Minimizer";
31} // namespace Prop
32
34namespace CostFunc {
35const static std::string UNWEIGHTED_LEAST_SQUARES = "Unweighted least squares";
36const static std::string WEIGHTED_LEAST_SQUARES = "Least squares";
37} // namespace CostFunc
38
40namespace Minimizer {
41const static std::string LEVENBERG_MARQUARDT_MD = "Levenberg-MarquardtMD";
42const static std::string LEVENBERG_MARQUARDT = "Levenberg-Marquardt";
43} // namespace Minimizer
44
51std::vector<double> filterRangesOutsideX(const std::vector<double> &ranges, const Mantid::API::MatrixWorkspace &ws,
52 const size_t wsIndex) {
53 const auto minX = ws.x(wsIndex).front();
54 const auto maxX = ws.x(wsIndex).back();
55 std::vector<double> filtered;
56 filtered.reserve(ranges.size());
57 for (size_t i = 0; i < ranges.size() / 2; ++i) {
58 const auto first = ranges[2 * i];
59 const auto second = ranges[2 * i + 1];
60 if (!(first > maxX || second < minX)) {
61 filtered.emplace_back(first);
62 filtered.emplace_back(second);
63 }
64 }
65 return filtered;
66}
67
74std::pair<double, double> totalRange(const std::vector<double> &ranges, const Mantid::API::MatrixWorkspace &ws,
75 const size_t wsIndex) {
76 const auto minX = ws.x(wsIndex).front();
77 const auto maxX = ws.x(wsIndex).back();
78 if (ranges.empty()) {
79 return std::pair<double, double>(minX, maxX);
80 }
81 const auto minmaxIt = std::minmax_element(ranges.cbegin(), ranges.cend());
82 const auto minEdge = *minmaxIt.first;
83 const auto maxEdge = *minmaxIt.second;
84 return std::pair<double, double>(std::min(minEdge, minX), std::max(maxEdge, maxX));
85}
86
92std::vector<double> includedRanges(const std::vector<double> &ranges, const std::pair<double, double> &totalRange) {
93 if (ranges.empty()) {
94 return {totalRange.first, totalRange.second};
95 }
96 // Sort the range edges keeping the information whether the edge
97 // 'starts' or 'ends' a range.
98 enum class Edge { start = -1, end = 1 };
99 std::vector<std::pair<double, Edge>> edges(ranges.size());
100 for (size_t i = 0; i < ranges.size(); ++i) {
101 edges[i].first = ranges[i];
102 edges[i].second = i % 2 == 0 ? Edge::start : Edge::end;
103 }
104 std::sort(edges.begin(), edges.end(), [](const std::pair<double, Edge> &p1, const std::pair<double, Edge> &p2) {
105 if (p1.first == p2.first)
106 return p1.second < p2.second;
107 return p1.first < p2.first;
108 });
109 // If an 'end' edge is followed by a 'start', we have a new range. Everything
110 // else can be merged.
111 std::vector<double> mergedRanges;
112 mergedRanges.reserve(ranges.size());
113 auto edgeIt = edges.begin();
114 mergedRanges.emplace_back(std::max(edges.front().first, totalRange.first));
115 while (edgeIt != edges.end()) {
116 auto endEdgeIt = edgeIt + 1;
117 while (endEdgeIt != edges.end()) {
118 const auto val = *endEdgeIt;
119 const auto prevVal = *(endEdgeIt - 1);
120 if (val.second == Edge::start && prevVal.second == Edge::end) {
121 mergedRanges.emplace_back(prevVal.first);
122 mergedRanges.emplace_back(val.first);
123 edgeIt = endEdgeIt;
124 break;
125 }
126 ++endEdgeIt;
127 }
128 ++edgeIt;
129 }
130 mergedRanges.emplace_back(std::min(edges.back().first, totalRange.second));
131 return mergedRanges;
132}
133
140std::vector<double> histogramRanges(const std::vector<double> &ranges, const Mantid::API::MatrixWorkspace &ws,
141 const size_t wsIndex) {
142 const auto filteredRanges = filterRangesOutsideX(ranges, ws, wsIndex);
143 if (!ranges.empty() && filteredRanges.empty()) {
144 throw std::runtime_error("The given XRanges mismatch with the histogram at workspace index " +
145 std::to_string(wsIndex));
146 }
147 const auto fullRange = totalRange(filteredRanges, ws, wsIndex);
148 return includedRanges(filteredRanges, fullRange);
149}
150
155std::vector<double> invertRanges(const std::vector<double> &ranges) {
156 std::vector<double> inversion(ranges.size() - 2);
157 for (size_t i = 1; i < ranges.size() - 1; ++i) {
158 inversion[i - 1] = ranges[i];
159 }
160 return inversion;
161}
162
173std::vector<double> executeFit(Mantid::API::Algorithm &fit, const std::string &function,
174 Mantid::API::MatrixWorkspace_sptr &ws, const size_t wsIndex,
175 const std::vector<double> &ranges, const std::string &costFunction,
176 const std::string &minimizer) {
177 const auto fitRanges = histogramRanges(ranges, *ws, wsIndex);
178 const auto excludedRanges = invertRanges(fitRanges);
179 fit.setProperty("Function", function);
180 fit.setProperty("InputWorkspace", ws);
181 fit.setProperty("WorkspaceIndex", static_cast<int>(wsIndex));
182 fit.setProperty("StartX", fitRanges.front());
183 fit.setProperty("EndX", fitRanges.back());
184 fit.setProperty("Exclude", excludedRanges);
185 fit.setProperty("Minimizer", minimizer);
186 fit.setProperty(Prop::COST_FUNCTION, costFunction);
187 fit.setProperty("CreateOutput", true);
188 fit.executeAsChildAlg();
189 Mantid::API::ITableWorkspace_sptr fitResult = fit.getProperty("OutputParameters");
190 std::vector<double> parameters(fitResult->rowCount() - 1);
191 for (size_t row = 0; row < parameters.size(); ++row) {
192 parameters[row] = fitResult->cell<double>(row, 1);
193 }
194 return parameters;
195}
196
202std::string makeFunctionString(const std::string &name, const std::vector<double> &parameters) {
203 const auto degree = parameters.size() - 1;
204 std::ostringstream s;
205 s << "name=" << name;
206 if (degree > 2)
207 s << ",n=" << degree;
208 for (size_t d = 0; d <= degree; ++d) {
209 s << ',' << 'A' << d << '=' << parameters[d];
210 }
211 return s.str();
212}
213
218std::string makeNameString(const size_t degree) {
219 std::ostringstream name;
220 switch (degree) {
221 case 0:
222 name << "FlatBackground";
223 break;
224 case 1:
225 name << "LinearBackground";
226 break;
227 case 2:
228 name << "Quadratic";
229 break;
230 default:
231 name << "Polynomial";
232 }
233 return name.str();
234}
235
242void evaluateInPlace(const std::string &name, const std::vector<double> &parameters, Mantid::API::MatrixWorkspace &ws,
243 const size_t wsIndex) {
244 const auto degree = parameters.size() - 1;
245 auto bkg = std::dynamic_pointer_cast<Mantid::API::IFunction1D>(
246 Mantid::API::FunctionFactory::Instance().createFunction(name));
247 if (degree > 2) {
248 Mantid::API::IFunction1D::Attribute att = bkg->getAttribute("n");
249 att.fromString(std::to_string(degree));
250 bkg->setAttribute("n", att);
251 }
252 for (size_t d = 0; d <= degree; ++d) {
253 std::string param = 'A' + std::to_string(d);
254 bkg->setParameter(param, parameters[d]);
255 }
256 auto *y = const_cast<double *>(ws.mutableY(wsIndex).rawData().data());
257 bkg->function1D(y, ws.points(wsIndex).rawData().data(), ws.y(wsIndex).size());
258}
259} // namespace
260
261namespace Mantid::Algorithms {
262
263// Register the algorithm into the AlgorithmFactory
264DECLARE_ALGORITHM(CalculatePolynomialBackground)
265
266//----------------------------------------------------------------------------------------------
267
268
269const std::string CalculatePolynomialBackground::name() const { return "CalculatePolynomialBackground"; }
270
273
276 return "CorrectionFunctions\\BackgroundCorrections";
277}
278
280const std::string CalculatePolynomialBackground::summary() const {
281 return "Fits a polynomial background to a workspace.";
282}
283
284//----------------------------------------------------------------------------------------------
288 auto increasingAxis = std::make_shared<API::IncreasingAxisValidator>();
289 auto nonnegativeInt = std::make_shared<Kernel::BoundedValidator<int>>();
290 nonnegativeInt->setLower(0);
291 auto orderedPairs = std::make_shared<Kernel::ArrayOrderedPairsValidator<double>>();
293 Prop::INPUT_WS, "", Kernel::Direction::Input, increasingAxis),
294 "An input workspace.");
296 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>(Prop::OUTPUT_WS, "", Kernel::Direction::Output),
297 "A workspace containing the fitted background.");
298 declareProperty(Prop::POLY_DEGREE, 0, nonnegativeInt, "Degree of the fitted polynomial.");
299 declareProperty(std::make_unique<Kernel::ArrayProperty<double>>(Prop::XRANGES, std::vector<double>(), orderedPairs),
300 "A list of fitting ranges given as pairs of X values.");
301 std::array<std::string, 2> costFuncOpts{{CostFunc::WEIGHTED_LEAST_SQUARES, CostFunc::UNWEIGHTED_LEAST_SQUARES}};
302 declareProperty(Prop::COST_FUNCTION, CostFunc::WEIGHTED_LEAST_SQUARES,
303 std::make_shared<Kernel::ListValidator<std::string>>(costFuncOpts),
304 "The cost function to be passed to the Fit algorithm.");
305 std::array<std::string, 2> minimizerOpts{{Minimizer::LEVENBERG_MARQUARDT_MD, Minimizer::LEVENBERG_MARQUARDT}};
306 declareProperty(Prop::MINIMIZER, Minimizer::LEVENBERG_MARQUARDT_MD,
307 std::make_shared<Kernel::ListValidator<std::string>>(minimizerOpts),
308 "The minimizer to be passed to the Fit algorithm.");
309}
310
311//----------------------------------------------------------------------------------------------
315 API::MatrixWorkspace_sptr inWS = getProperty(Prop::INPUT_WS);
316
317 API::MatrixWorkspace_sptr outWS{DataObjects::create<DataObjects::Workspace2D>(*inWS)};
318 const std::vector<double> inputRanges = getProperty(Prop::XRANGES);
319 const std::string costFunction = getProperty(Prop::COST_FUNCTION);
320 const std::string minimizer = getProperty(Prop::MINIMIZER);
321 const auto polyDegree = static_cast<size_t>(static_cast<int>(getProperty(Prop::POLY_DEGREE)));
322 const std::vector<double> initialParams(polyDegree + 1, 0.1);
323 const auto name = makeNameString(polyDegree);
324 const auto fitFunction = makeFunctionString(name, initialParams);
325 const auto nHistograms = static_cast<int64_t>(inWS->getNumberHistograms());
326 API::Progress progress(this, 0, 1.0, nHistograms);
327 PARALLEL_FOR_IF(Kernel::threadSafe(*inWS, *outWS))
328 for (int64_t i = 0; i < nHistograms; ++i) {
330 const bool logging{false};
331 auto fit = createChildAlgorithm("Fit", 0, 0, logging);
332 const auto parameters = executeFit(*fit, fitFunction, inWS, i, inputRanges, costFunction, minimizer);
333 evaluateInPlace(name, parameters, *outWS, i);
334 progress.report();
336 }
338
339 setProperty(Prop::OUTPUT_WS, outWS);
340}
341
342} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
void executeAsChildAlg() override
Execute as a Child Algorithm.
Definition: Algorithm.cpp:769
Attribute is a non-fitting parameter.
Definition: IFunction.h:282
void fromString(const std::string &str)
Set value from a string.
Definition: IFunction.cpp:971
Base MatrixWorkspace Abstract Class.
HistogramData::Points points(const size_t index) const
const HistogramData::HistogramX & x(const size_t index) const
HistogramData::HistogramY & mutableY(const size_t index) &
const HistogramData::HistogramY & y(const size_t index) const
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
CalculatePolynomialBackground : This algorithm fits a polynomial background to an input workspace and...
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
int version() const override
Algorithm's version for identification.
const std::string category() const override
Algorithm's category for identification.
void init() override
Initialize the algorithm's properties.
const std::string name() const override
Algorithms name for identification.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
ListValidator is a validator that requires the value of a property to be one of a defined list of pos...
Definition: ListValidator.h:29
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
String constants for cost function options.
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
boost::graph_traits< PeakGraph >::edge_descriptor Edge
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
String constants for minimizer options.
String constants for algorithm's properties.
STL namespace.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54