Mantid
Loading...
Searching...
No Matches
Rebin.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#include "MantidHistogramData/Exception.h"
9#include "MantidHistogramData/Rebin.h"
10
11#include "MantidAPI/Axis.h"
21
22namespace Mantid::Algorithms {
23
24// Register the class into the algorithm factory
26
27using namespace Kernel;
28using namespace API;
29using DataObjects::EventList;
30using DataObjects::EventWorkspace;
33using HistogramData::BinEdges;
34using HistogramData::Frequencies;
35using HistogramData::FrequencyStandardDeviations;
36using HistogramData::Histogram;
37using HistogramData::Exception::InvalidBinEdgesError;
38
39//---------------------------------------------------------------------------------------------
40// Public static methods
41//---------------------------------------------------------------------------------------------
42
50std::vector<double> Rebin::rebinParamsFromInput(const std::vector<double> &inParams,
51 const API::MatrixWorkspace &inputWS, Kernel::Logger &logger) {
52 std::vector<double> rbParams;
53 // The validator only passes parameters with size 1, or 3xn. No need to check again here
54 if (inParams.size() >= 3) {
55 // Input are min, delta, max
56 rbParams = inParams;
57 } else if (inParams.size() == 1) {
58 double xmin = 0.;
59 double xmax = 0.;
60 inputWS.getXMinMax(xmin, xmax);
61
62 logger.information() << "Using the current min and max as default " << xmin << ", " << xmax << '\n';
63 rbParams.resize(3);
64 rbParams[0] = xmin;
65 rbParams[1] = inParams[0];
66 rbParams[2] = xmax;
67 if ((rbParams[1] < 0.) && (xmin < 0.) && (xmax > 0.)) {
68 std::stringstream msg;
69 msg << "Cannot create logarithmic binning that changes sign (xmin=" << xmin << ", xmax=" << xmax << ")";
70 throw std::runtime_error(msg.str());
71 }
72 }
73 return rbParams;
74}
75
76//---------------------------------------------------------------------------------------------
77// Public methods
78//---------------------------------------------------------------------------------------------
79
81std::map<std::string, std::string> Rebin::validateInputs() {
82 std::map<std::string, std::string> helpMessages;
83 if (existsProperty("Power") && !isDefault("Power")) {
84 const double power = getProperty("Power");
85
86 // attempt to roughly guess how many bins these parameters imply
87 double roughEstimate = 0;
88
89 if (!isDefault("Params")) {
90 const std::vector<double> params = getProperty("Params");
91
92 // Five significant places of the Euler-Mascheroni constant is probably more than enough for our needs
93 double eulerMascheroni = 0.57721;
94
95 // Params is check by the validator first, so we can assume it is in a correct format
96 for (size_t i = 0; i < params.size() - 2; i += 2) {
97 double upperLimit = params[i + 2];
98 double lowerLimit = params[i];
99 double factor = params[i + 1];
100
101 if (factor <= 0) {
102 helpMessages["Params"] = "Provided width value cannot be negative for inverse power binning.";
103 return helpMessages;
104 }
105
106 if (power == 1) {
107 roughEstimate += std::exp((upperLimit - lowerLimit) / factor - eulerMascheroni);
108 } else {
109 roughEstimate += std::pow(((upperLimit - lowerLimit) / factor) * (1 - power) + 1, 1 / (1 - power));
110 }
111 }
112 }
113
114 // Prevent the user form creating too many bins
115 if (roughEstimate > 10000) {
116 helpMessages["Power"] = "This binning is expected to give more than 10000 bins.";
117 }
118 }
119 return helpMessages;
120}
121
126 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
127 "Workspace containing the input data");
128 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
129 "The name to give the output workspace");
130
131 declareProperty(std::make_unique<ArrayProperty<double>>("Params", std::make_shared<RebinParamsValidator>()),
132 "A comma separated list of first bin boundary, width, last bin boundary. "
133 "Optionally this can be followed by a comma and more widths and last boundary pairs. "
134 "Optionally this can also be a single number, which is the bin width. In this case, the boundary of "
135 "binning will be determined by minimum and maximum TOF values among all events, or previous binning "
136 "boundary, in case of event Workspace, or non-event Workspace, respectively. "
137 "Negative width values indicate logarithmic binning.");
138
139 declareProperty("PreserveEvents", true,
140 "Keep the output workspace as an EventWorkspace, if the input has events. If the input and output "
141 "EventWorkspace names are the same, only the X bins are set, which is very quick. If false, then the "
142 "workspace gets converted to a Workspace2D histogram.");
143
144 declareProperty("FullBinsOnly", false, "Omit the final bin if its width is smaller than the step size");
145
146 declareProperty("IgnoreBinErrors", false,
147 "Ignore errors related to zero/negative bin widths in input/output workspaces. When ignored, the "
148 "signal and errors are set to zero");
149
151 "UseReverseLogarithmic", false,
152 "For logarithmic intervals, the splitting starts from the end and goes back to the start, ie the bins are bigger "
153 "at the start getting exponentially smaller until they reach the end. For these bins, the FullBinsOnly flag is "
154 "ignored.");
155
156 auto powerValidator = std::make_shared<Mantid::Kernel::BoundedValidator<double>>();
157 powerValidator->setLower(0);
158 powerValidator->setUpper(1);
159 declareProperty("Power", 0., powerValidator,
160 "Splits the interval in bins which actual width is equal to requested width / (i ^ power); default "
161 "is linear. Power must be between 0 and 1.");
162}
163
170 // Get the input workspace
171 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
172 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
173
174 // Are we preserving event workspace-iness?
175 bool PreserveEvents = getProperty("PreserveEvents");
176
177 // Rebinning in-place
178 bool inPlace = (inputWS == outputWS);
179
180 std::vector<double> rbParams = rebinParamsFromInput(getProperty("Params"), *inputWS, g_log);
181
182 const bool dist = inputWS->isDistribution();
183 const bool isHist = inputWS->isHistogramData();
184
185 // workspace independent determination of length
186 const auto histnumber = static_cast<int>(inputWS->getNumberHistograms());
187
188 bool fullBinsOnly = getProperty("FullBinsOnly");
189 bool useReverseLog = getProperty("UseReverseLogarithmic");
190 double power = getProperty("Power");
191
192 double xmin = 0.;
193 double xmax = 0.;
194 inputWS->getXMinMax(xmin, xmax);
195
196 HistogramData::BinEdges XValues_new(0);
197 // create new output X axis
198 static_cast<void>(VectorHelper::createAxisFromRebinParams(rbParams, XValues_new.mutableRawData(), true, fullBinsOnly,
199 xmin, xmax, useReverseLog, power));
200
201 // Now, determine if the input workspace is actually an EventWorkspace
202 EventWorkspace_const_sptr eventInputWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWS);
203
204 if (eventInputWS != nullptr) {
205 //------- EventWorkspace as input -------------------------------------
206
207 if (PreserveEvents) {
208 if (!inPlace) {
209 outputWS = inputWS->clone();
210 }
211 auto eventOutputWS = std::dynamic_pointer_cast<EventWorkspace>(outputWS);
212 // This only sets the X axis. Actual rebinning will be done upon data
213 // access.
214 eventOutputWS->setAllX(XValues_new);
215 } else {
216 //--------- Different output, OR you're inplace but not preserving Events
217 g_log.information() << "Creating a Workspace2D from the EventWorkspace " << eventInputWS->getName() << ".\n";
218 outputWS = DataObjects::create<DataObjects::Workspace2D>(*inputWS, histnumber, XValues_new);
219
220 // Initialize progress reporting.
221 Progress prog(this, 0.0, 1.0, histnumber);
222
223 // Go through all the histograms and set the data
224 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
225 for (int i = 0; i < histnumber; ++i) {
227 // Get a const event list reference. eventInputWS->dataY() doesn't work.
228 const EventList &el = eventInputWS->getSpectrum(i);
229 MantidVec y_data, e_data;
230 // The EventList takes care of histogramming.
231 el.generateHistogram(XValues_new.rawData(), y_data, e_data);
232
233 // Copy the data over.
234 outputWS->mutableY(i) = y_data;
235 outputWS->mutableE(i) = e_data;
236
237 // Report progress
238 prog.report(name());
240 }
242
243 // Copy all the axes
244 for (int i = 1; i < inputWS->axes(); i++) {
245 outputWS->replaceAxis(i, std::unique_ptr<Axis>(inputWS->getAxis(i)->clone(outputWS.get())));
246 outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit();
247 }
248
249 // Copy the units over too.
250 for (int i = 0; i < outputWS->axes(); ++i)
251 outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit();
252 outputWS->setYUnit(eventInputWS->YUnit());
253 outputWS->setYUnitLabel(eventInputWS->YUnitLabel());
254 }
255
256 // Assign it to the output workspace property
257 setProperty("OutputWorkspace", outputWS);
258
259 } // END ---- EventWorkspace
260
261 else
262
263 { //------- Workspace2D or other MatrixWorkspace ---------------------------
264
265 if (!isHist) {
266 g_log.information() << "Rebin: Converting Data to Histogram.\n";
267 Mantid::API::Algorithm_sptr ChildAlg = createChildAlgorithm("ConvertToHistogram");
268 ChildAlg->initialize();
269 ChildAlg->setProperty("InputWorkspace", inputWS);
270 ChildAlg->execute();
271 inputWS = ChildAlg->getProperty("OutputWorkspace");
272 }
273
274 // make output Workspace the same type is the input, but with new length of
275 // signal array
276 outputWS = DataObjects::create<API::HistoWorkspace>(*inputWS, histnumber, XValues_new);
277
278 // Copy over the 'vertical' axis
279 if (inputWS->axes() > 1)
280 outputWS->replaceAxis(1, std::unique_ptr<Axis>(inputWS->getAxis(1)->clone(outputWS.get())));
281 bool ignoreBinErrors = getProperty("IgnoreBinErrors");
282
283 Progress prog(this, 0.0, 1.0, histnumber);
284 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
285 for (int hist = 0; hist < histnumber; ++hist) {
287
288 try {
289 outputWS->setHistogram(hist, HistogramData::rebin(inputWS->histogram(hist), XValues_new));
290 } catch (InvalidBinEdgesError &) {
291 if (ignoreBinErrors)
292 outputWS->setBinEdges(hist, XValues_new);
293 else
294 throw;
295 }
296 prog.report(name());
298 }
300 outputWS->setDistribution(dist);
301
302 // Now propagate any masking correctly to the output workspace
303 // More efficient to have this in a separate loop because
304 // MatrixWorkspace::maskBins blocks multi-threading
305 for (int i = 0; i < histnumber; ++i) {
306 if (inputWS->hasMaskedBins(i)) // Does the current spectrum have any masked bins?
307 {
308 this->propagateMasks(inputWS, outputWS, i);
309 }
310 }
311 // Copy the units over too.
312 for (int i = 0; i < outputWS->axes(); ++i) {
313 outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit();
314 }
315
316 if (!isHist) {
317 g_log.information() << "Rebin: Converting Data back to Data Points.\n";
318 Mantid::API::Algorithm_sptr ChildAlg = createChildAlgorithm("ConvertToPointData");
319 ChildAlg->initialize();
320 ChildAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", outputWS);
321 ChildAlg->execute();
322 outputWS = ChildAlg->getProperty("OutputWorkspace");
323 }
324
325 // Assign it to the output workspace property
326 setProperty("OutputWorkspace", outputWS);
327
328 } // END ---- Workspace2D
329}
330
340 int hist) {
341 // Not too happy with the efficiency of this way of doing it, but it's a lot
342 // simpler to use the
343 // existing rebin algorithm to distribute the weights than to re-implement it
344 // for this
345
346 MantidVec masked_bins, weights;
347 // Get a reference to the list of masked bins for this spectrum
348 const MatrixWorkspace::MaskList &mask = inputWS->maskedBins(hist);
349 // Now iterate over the list, building up a vector of the masked bins
350 auto it = mask.cbegin();
351 auto &XValues = inputWS->x(hist);
352 masked_bins.emplace_back(XValues[(*it).first]);
353 weights.emplace_back((*it).second);
354 masked_bins.emplace_back(XValues[(*it).first + 1]);
355 for (++it; it != mask.end(); ++it) {
356 const double currentX = XValues[(*it).first];
357 // Add an intermediate bin with zero weight if masked bins aren't
358 // consecutive
359 if (masked_bins.back() != currentX) {
360 weights.emplace_back(0.0);
361 masked_bins.emplace_back(currentX);
362 }
363 weights.emplace_back((*it).second);
364 masked_bins.emplace_back(XValues[(*it).first + 1]);
365 }
366
368 auto errSize = weights.size();
369 Histogram oldHist(BinEdges(std::move(masked_bins)), Frequencies(std::move(weights)),
370 FrequencyStandardDeviations(errSize, 0));
371 // Use rebin function to redistribute the weights. Note that distribution flag
372 // is set
373 bool ignoreErrors = getProperty("IgnoreBinErrors");
374
375 try {
376 auto newHist = HistogramData::rebin(oldHist, outputWS->binEdges(hist));
377 auto &newWeights = newHist.y();
378
379 // Now process the output vector and fill the new masking list
380 for (size_t index = 0; index < newWeights.size(); ++index) {
381 if (newWeights[index] > 0.0)
382 outputWS->flagMasked(hist, index, newWeights[index]);
383 }
384 } catch (InvalidBinEdgesError &) {
385 if (!ignoreErrors)
386 throw;
387 }
388}
389
390} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#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.
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
bool existsProperty(const std::string &name) const override
Checks whether the named property is already in the list of managed property.
Definition: Algorithm.cpp:2008
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
bool isDefault(const std::string &name) const
Definition: Algorithm.cpp:2084
Base MatrixWorkspace Abstract Class.
std::map< size_t, double > MaskList
Masked bins for each spectrum are stored as a set of pairs containing <bin index, weight>
virtual void getXMinMax(double &xmin, double &xmax) const
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void propagateMasks(const API::MatrixWorkspace_const_sptr &inputWS, const API::MatrixWorkspace_sptr &outputWS, int hist)
Takes the masks in the input workspace and apportions the weights into the new bins that overlap with...
Definition: Rebin.cpp:339
void exec() override
Executes the rebin algorithm.
Definition: Rebin.cpp:169
void init() override
Initialisation method.
Definition: Rebin.cpp:125
const std::string name() const override
Algorithm's name for identification overriding a virtual method.
Definition: Rebin.h:41
static std::vector< double > rebinParamsFromInput(const std::vector< double > &inParams, const API::MatrixWorkspace &inputWS, Kernel::Logger &logger)
Return the rebin parameters from a user input.
Definition: Rebin.cpp:50
std::map< std::string, std::string > validateInputs() override
Validate that the input properties are sane.
Definition: Rebin.cpp:81
A class for holding :
Definition: EventList.h:56
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.
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< Algorithm > Algorithm_sptr
Typedef for a shared pointer to an Algorithm.
Definition: Algorithm.h:61
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
int MANTID_KERNEL_DLL createAxisFromRebinParams(const std::vector< double > &params, std::vector< double > &xnew, const bool resize_xnew=true, const bool full_bins_only=false, const double xMinHint=std::nan(""), const double xMaxHint=std::nan(""), const bool useReverseLogarithmic=false, const double power=-1)
Creates a new output X array given a 'standard' set of rebinning parameters.
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
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54