/
MaskNonOverlappingBins.cpp
234 lines (214 loc) · 9.95 KB
/
MaskNonOverlappingBins.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAlgorithms/MaskNonOverlappingBins.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/Unit.h"
namespace {
/// Constants for the algorithm's property names.
namespace Prop {
std::string const CHECK_SORTING{"CheckSortedX"};
std::string const COMPARISON_WS{"ComparisonWorkspace"};
std::string const INPUT_WS{"InputWorkspace"};
std::string const MASK_PARTIAL{"MaskPartiallyOverlapping"};
std::string const OUTPUT_WS{"OutputWorkspace"};
std::string const RAGGEDNESS{"RaggedInputs"};
} // namespace Prop
/// Constants for the RaggedInputs property.
namespace Raggedness {
std::string const CHECK{"Check"};
std::string const RAGGED{"Ragged"};
std::string const NONRAGGED{"Common Bins"};
} // namespace Raggedness
/// Return true if X data is sorted in ascending order.
bool isXSorted(Mantid::API::MatrixWorkspace const &ws) {
int unsorted{0};
PARALLEL_FOR_IF(Mantid::Kernel::threadSafe(ws))
for (int64_t i = 0; i < static_cast<int64_t>(ws.getNumberHistograms()); ++i) {
auto const &Xs = ws.x(i);
if (!std::is_sorted(Xs.cbegin(), Xs.cend())) {
PARALLEL_ATOMIC
++unsorted;
}
}
return unsorted == 0;
}
/// Holds limiting bin indices for masking
struct BinIndices {
size_t frontEndIndex; ///< Mask from 0 to this index.
size_t backBeginIndex; ///< Mask from this index to size.
};
/// Return the proper masking limits for non-overlapping bins.
BinIndices maskingLimits(Mantid::API::MatrixWorkspace const &ws, Mantid::API::MatrixWorkspace const &comparisonWS,
bool const maskPartial, size_t const histogramIndex) {
auto const &Xs = ws.x(histogramIndex);
auto const &comparisonXs = comparisonWS.x(histogramIndex);
// At the moment we only support increasing X.
auto const startX = comparisonXs.front();
auto const endX = comparisonXs.back();
// There is no Y corresponding to the last bin edge,
// therefore iterate only to cend() - 1.
auto frontEnd = std::lower_bound(Xs.cbegin(), Xs.cend() - 1, startX);
if (!maskPartial && frontEnd != Xs.cbegin() && *frontEnd != startX) {
--frontEnd;
}
auto backBegin = std::lower_bound(frontEnd, Xs.cend() - 1, endX);
if (maskPartial && backBegin != Xs.cbegin() && *backBegin > endX) {
--backBegin;
}
BinIndices limits;
limits.frontEndIndex = static_cast<size_t>(std::distance(Xs.cbegin(), frontEnd));
limits.backBeginIndex = static_cast<size_t>(std::distance(Xs.cbegin(), backBegin));
return limits;
}
/// Mask the start and end of histogram at histogram index.
void maskBinsWithinLimits(Mantid::API::MatrixWorkspace &ws, size_t const histogramIndex, BinIndices const &limits) {
auto const xSize = ws.x(histogramIndex).size();
for (size_t binIndex = 0; binIndex < limits.frontEndIndex; ++binIndex) {
ws.flagMasked(histogramIndex, binIndex);
}
for (size_t binIndex = limits.backBeginIndex; binIndex < xSize - 1; ++binIndex) {
ws.flagMasked(histogramIndex, binIndex);
}
}
} // namespace
namespace Mantid::Algorithms {
DECLARE_ALGORITHM(MaskNonOverlappingBins)
/// Algorithms name for identification. @see Algorithm::name
std::string const MaskNonOverlappingBins::name() const { return "MaskNonOverlappingBins"; }
/// Algorithm's version for identification. @see Algorithm::version
int MaskNonOverlappingBins::version() const { return 1; }
/// Algorithm's category for identification. @see Algorithm::category
std::string const MaskNonOverlappingBins::category() const { return "Transforms\\Masking"; }
/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
std::string const MaskNonOverlappingBins::summary() const {
return "Marks bins in InputWorkspace which are out of the X range of the "
"second workspace.";
}
/// Return a list of the names of related algorithms.
std::vector<std::string> const MaskNonOverlappingBins::seeAlso() const { return {"MaskBins", "MaskBinsIf"}; }
/** Initialize the algorithm's properties.
*/
void MaskNonOverlappingBins::init() {
declareProperty(std::make_unique<API::WorkspaceProperty<>>(Prop::INPUT_WS, "", Kernel::Direction::Input),
"A workspace to mask.");
declareProperty(std::make_unique<API::WorkspaceProperty<>>(Prop::OUTPUT_WS, "", Kernel::Direction::Output),
"The masked workspace.");
declareProperty(std::make_unique<API::WorkspaceProperty<>>(Prop::COMPARISON_WS, "", Kernel::Direction::Input),
"A workspace to compare the InputWorkspace's binning to.");
declareProperty(Prop::MASK_PARTIAL, false, "If true, mask also bins that overlap only partially.");
std::vector<std::string> const options{Raggedness::CHECK, Raggedness::RAGGED, Raggedness::NONRAGGED};
auto raggednessOptions = std::make_shared<Kernel::ListValidator<std::string>>(options);
declareProperty(Prop::RAGGEDNESS, Raggedness::CHECK, raggednessOptions,
"Choose whether the input workspaces have common bins, are "
"ragged, or if the algorithm should check.");
declareProperty(Prop::CHECK_SORTING, true,
"If true, the algorithm ensures that both workspaces have X sorted in "
"ascending order.");
}
/// Returns a map from property name to message in case of invalid values.
std::map<std::string, std::string> MaskNonOverlappingBins::validateInputs() {
std::map<std::string, std::string> issues;
API::MatrixWorkspace_const_sptr inputWS = getProperty(Prop::INPUT_WS);
API::MatrixWorkspace_const_sptr comparisonWS = getProperty(Prop::COMPARISON_WS);
if (!inputWS) {
issues[Prop::INPUT_WS] = "The " + Prop::INPUT_WS + " must be a MatrixWorkspace.";
}
if (!comparisonWS) {
issues[Prop::COMPARISON_WS] = "The " + Prop::COMPARISON_WS + " must be a MatrixWorkspace.";
}
if (!issues.empty()) {
return issues;
}
if (inputWS->getNumberHistograms() != comparisonWS->getNumberHistograms()) {
issues[Prop::COMPARISON_WS] = "The number of histogams mismatches with " + Prop::INPUT_WS;
}
if (!inputWS->isHistogramData()) {
issues[Prop::INPUT_WS] = "The workspace contains point data, not histograms.";
}
if (!comparisonWS->isHistogramData()) {
issues[Prop::COMPARISON_WS] = "The workspace contains point data, not histograms.";
}
auto const inputAxis = inputWS->getAxis(0);
auto const comparisonAxis = comparisonWS->getAxis(0);
if (inputAxis && comparisonAxis) {
if (*(inputAxis->unit()) != *(comparisonAxis->unit())) {
issues[Prop::COMPARISON_WS] = "X units do not match with " + Prop::INPUT_WS;
}
}
return issues;
}
/** Execute the algorithm.
*/
void MaskNonOverlappingBins::exec() {
API::MatrixWorkspace_sptr inputWS = getProperty(Prop::INPUT_WS);
API::MatrixWorkspace_sptr outputWS = getProperty(Prop::OUTPUT_WS);
if (inputWS != outputWS) {
outputWS = inputWS->clone();
}
API::MatrixWorkspace_const_sptr comparisonWS = getProperty(Prop::COMPARISON_WS);
checkXSorting(*inputWS, *comparisonWS);
if (isCommonBins(*inputWS, *comparisonWS)) {
processNonRagged(*inputWS, *comparisonWS, *outputWS);
} else {
processRagged(*inputWS, *comparisonWS, *outputWS);
}
setProperty(Prop::OUTPUT_WS, outputWS);
}
/// Throw if the workspaces don't have sorted X.
void MaskNonOverlappingBins::checkXSorting(API::MatrixWorkspace const &inputWS,
API::MatrixWorkspace const &comparisonWS) {
bool const checkSorting = getProperty(Prop::CHECK_SORTING);
if (checkSorting) {
if (!isXSorted(inputWS)) {
throw std::invalid_argument(Prop::INPUT_WS + " has unsorted X.");
}
if (!isXSorted(comparisonWS)) {
throw std::invalid_argument(Prop::COMPARISON_WS + " has unsorted X.");
}
}
}
/// Return true if the workspaces should be considered as having common bins.
bool MaskNonOverlappingBins::isCommonBins(API::MatrixWorkspace const &inputWS,
API::MatrixWorkspace const &comparisonWS) {
std::string const choice = getProperty(Prop::RAGGEDNESS);
if (choice == Raggedness::CHECK) {
return inputWS.isCommonBins() && comparisonWS.isCommonBins();
} else {
return choice == Raggedness::NONRAGGED;
}
}
/// Mask all types of workspaces, ragged or nonragged.
void MaskNonOverlappingBins::processRagged(API::MatrixWorkspace const &inputWS,
API::MatrixWorkspace const &comparisonWS, API::MatrixWorkspace &outputWS) {
bool const maskPartial = getProperty(Prop::MASK_PARTIAL);
auto const nHist = inputWS.getNumberHistograms();
API::Progress progress(this, 0., 1., nHist);
// Tried to parallelize this loop but the performance test showed
// regression. Hence no PARALLEL_FOR_IF.
for (size_t histogramIndex = 0; histogramIndex < nHist; ++histogramIndex) {
auto const limits = maskingLimits(inputWS, comparisonWS, maskPartial, histogramIndex);
maskBinsWithinLimits(outputWS, histogramIndex, limits);
progress.report("Masking nonoverlapping bins");
}
}
/// Mask only workspace with same X in all histograms.
void MaskNonOverlappingBins::processNonRagged(API::MatrixWorkspace const &inputWS,
API::MatrixWorkspace const &comparisonWS,
API::MatrixWorkspace &outputWS) {
bool const maskPartial = getProperty(Prop::MASK_PARTIAL);
auto const nHist = inputWS.getNumberHistograms();
API::Progress progress(this, 0., 1., nHist);
auto const limits = maskingLimits(inputWS, comparisonWS, maskPartial, 0);
for (size_t histogramIndex = 0; histogramIndex < nHist; ++histogramIndex) {
maskBinsWithinLimits(outputWS, histogramIndex, limits);
progress.report("Masking nonoverlapping bins");
}
}
} // namespace Mantid::Algorithms