Mantid
Loading...
Searching...
No Matches
ResampleX.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 "MantidAPI/Axis.h"
16
17#include <sstream>
18
19namespace Mantid::Algorithms {
20using namespace API;
21using namespace DataObjects;
22using namespace Kernel;
23using HistogramData::BinEdges;
24using std::map;
25using std::string;
26using std::stringstream;
27using std::vector;
28
29// Register the algorithm into the AlgorithmFactory
30DECLARE_ALGORITHM(ResampleX)
31
32//----------------------------------------------------------------------------------------------
34const std::string ResampleX::name() const { return "ResampleX"; }
35
37int ResampleX::version() const { return 1; }
38
40const std::string ResampleX::category() const { return "Transforms\\Rebin"; }
41
42const std::string ResampleX::alias() const { return ""; }
43
44//----------------------------------------------------------------------------------------------
48 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input), "An input workspace.");
49 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
50 "An output workspace.");
51
52 declareProperty(std::make_unique<ArrayProperty<double>>("XMin"),
53 "A comma separated list of the XMin for every spectrum. (Optional)");
54 declareProperty(std::make_unique<ArrayProperty<double>>("XMax"),
55 "A comma separated list of the XMax for every spectrum. (Optional)");
56
57 auto min = std::make_shared<BoundedValidator<int>>();
58 min->setLower(1);
59 declareProperty("NumberBins", 0, min, "Number of bins to split up each spectrum into.");
60 declareProperty("LogBinning", false, "Use logarithmic binning. If false use constant step sizes.");
61
62 declareProperty("PreserveEvents", true,
63 "Keep the output workspace as an EventWorkspace, if the "
64 "input has events (default).\n"
65 "If the input and output EventWorkspace names are the same, "
66 "only the X bins are set, which is very quick.\n"
67 "If false, then the workspace gets converted to a "
68 "Workspace2D histogram.");
69}
70
74map<string, string> ResampleX::validateInputs() {
75 map<string, string> errors;
76 vector<double> xmins = getProperty("XMin");
77 vector<double> xmaxs = getProperty("XMax");
78 if ((!xmins.empty()) && (!xmaxs.empty())) {
79 if (xmins.size() != xmaxs.size()) {
80 stringstream msg;
81 msg << "XMin and XMax do not define same number of spectra (" << xmins.size() << " != " << xmaxs.size() << ")";
82 errors.emplace("XMax", msg.str());
83 } else {
84 size_t size = xmins.size();
85 for (size_t i = 0; i < size; ++i) {
86 if (xmins[i] >= xmaxs[i]) {
87 stringstream msg;
88 msg << "XMin (" << xmins[i] << ") cannot be greater than XMax (" << xmaxs[i] << ")";
89 errors.emplace("XMax", msg.str());
90 }
91 }
92 }
93 }
94
95 return errors;
96}
97
110string determineXMinMax(const MatrixWorkspace_sptr &inputWS, vector<double> &xmins, vector<double> &xmaxs) {
111 const size_t numSpectra = inputWS->getNumberHistograms();
112
113 // pad out the ranges by copying the first value to the rest that are needed
114 if (xmins.size() == 1 && numSpectra > xmins.size()) {
115 const double value = xmins.front();
116 xmins.insert(xmins.end(), numSpectra - xmins.size(), value);
117 }
118 if (xmaxs.size() == 1 && numSpectra > xmaxs.size()) {
119 const double value = xmaxs.front();
120 xmaxs.insert(xmaxs.end(), numSpectra - xmaxs.size(), value);
121 }
122
123 // should the individiual values be calculated?
124 const bool updateXMins = xmins.empty(); // they weren't set
125 const bool updateXMaxs = xmaxs.empty(); // they weren't set
126
127 stringstream msg;
128
129 // determine overall xmin/xmax
130 double xmin_wksp = inputWS->getXMin();
131 double xmax_wksp = inputWS->getXMax();
132 EventWorkspace_const_sptr inputEventWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWS);
133 if (inputEventWS != nullptr && inputEventWS->getNumberEvents() > 0) {
134 xmin_wksp = inputEventWS->getTofMin();
135 xmax_wksp = inputEventWS->getTofMax();
136 }
137
138 for (size_t i = 0; i < numSpectra; ++i) {
139 // determine ranges if necessary
140 if (updateXMins || updateXMaxs) {
141 const auto &xvalues = inputWS->x(i);
142 if (updateXMins) {
143 const auto minimum = xvalues.front();
144 if (std::isnan(minimum) || minimum >= xmax_wksp) {
145 xmins.emplace_back(xmin_wksp);
146 } else {
147 xmins.emplace_back(minimum);
148 }
149 }
150 if (updateXMaxs) {
151 const auto maximum = xvalues.back();
152 if (std::isnan(maximum) || maximum <= xmin_wksp) {
153 xmaxs.emplace_back(xmax_wksp);
154 } else {
155 xmaxs.emplace_back(maximum);
156 }
157 }
158 }
159
160 // error check the ranges
161 if (xmins[i] >= xmaxs[i]) {
162 if (!msg.str().empty())
163 msg << ", ";
164 msg << "at wksp_index=" << i << " XMin >= XMax (" << xmins[i] << " >= " << xmaxs[i] << ")";
165 }
166 }
167
168 return msg.str(); // empty string means nothing went wrong
169}
170
179void ResampleX::setOptions(const int numBins, const bool useLogBins, const bool isDist) {
180 m_numBins = numBins;
181 m_useLogBinning = useLogBins;
182 m_isDistribution = isDist;
183}
184
194double ResampleX::determineBinning(MantidVec &xValues, const double xmin, const double xmax) {
195 xValues.clear(); // clear out the x-values
196
197 int numBoundaries(0);
198 int reqNumBoundaries(m_numBins);
199 int expNumBoundaries(m_numBins);
201 reqNumBoundaries -= 1; // to get the VectorHelper to do the right thing
202 else
203 expNumBoundaries += 1; // should be one more bin boundary for histograms
204
205 vector<double> params; // xmin, delta, xmax
206 params.emplace_back(xmin);
207 params.emplace_back(0.); // dummy delta value
208 params.emplace_back(xmax);
209
210 // constant binning is easy
211 if (m_useLogBinning) {
212 if (xmin == 0)
213 throw std::invalid_argument("Cannot calculate log of xmin=0");
214 if (xmax == 0)
215 throw std::invalid_argument("Cannot calculate log of xmax=0");
216 if (xmin < 0. && xmax > 0.) {
217 std::stringstream msg;
218 msg << "Cannot calculate logorithmic binning that changes sign (xmin=" << xmin << ", xmax=" << xmax << ")";
219 throw std::invalid_argument(msg.str());
220 }
221
222 const int MAX_ITER(100); // things went wrong if we get this far
223
224 // starting delta value assuming everything happens exactly
225 double delta = (log(xmax) - log(xmin)) / static_cast<double>(m_numBins);
226 double shift = .1;
227 int sign = 0;
228 for (int numIter = 0; numIter < MAX_ITER; ++numIter) {
229 params[1] = -1. * delta;
230 if (!m_isDistribution)
231 params[2] = xmax + delta;
232 numBoundaries = VectorHelper::createAxisFromRebinParams(params, xValues, true);
233
234 if (numBoundaries == expNumBoundaries) {
235 double diff = (xmax - xValues.back());
236 if (diff != 0.) {
237 g_log.debug() << "Didn't get the exact xmax value: [xmax - xValues.back()=" << diff
238 << "] [relative diff = " << fabs(100. * diff / xmax) << "%]\n";
239 g_log.debug() << "Resetting final x-value to xmax\n";
240 *(xValues.rbegin()) = xmax;
241 }
242 break;
243 } else if (numBoundaries > expNumBoundaries) // too few points
244 {
245 delta *= (1. + shift);
246 if (sign < 0)
247 shift *= .9;
248 sign = 1;
249 } else // too many points
250 {
251 delta *= (1. - shift);
252 if (sign > 0)
253 shift *= .9;
254 sign = -1;
255 }
256 }
257 } else {
258 params[1] = (xmax - xmin) / static_cast<double>(reqNumBoundaries);
259 numBoundaries = VectorHelper::createAxisFromRebinParams(params, xValues, true);
260 }
261
262 if (numBoundaries != expNumBoundaries) {
263 g_log.warning() << "Did not generate the requested number of bins: generated " << numBoundaries << " requested "
264 << expNumBoundaries << "(xmin=" << xmin << ", xmax=" << xmax << ")\n";
265 }
266
267 // return the delta value so the caller can do debug printing
268 return params[1];
269}
270
271//----------------------------------------------------------------------------------------------
275 // generically having access to the input workspace is a good idea
276 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
277 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
278 bool inPlace = (inputWS == outputWS); // Rebinning in-place
279 m_isDistribution = inputWS->isDistribution();
280 m_isHistogram = inputWS->isHistogramData();
281 const auto numSpectra = static_cast<int>(inputWS->getNumberHistograms());
282
283 // the easy parameters
284 m_useLogBinning = getProperty("LogBinning");
285 m_numBins = getProperty("NumberBins");
286 m_preserveEvents = getProperty("PreserveEvents");
287
288 // determine the xmin/xmax for the workspace
289 vector<double> xmins = getProperty("XMin");
290 vector<double> xmaxs = getProperty("XMax");
291 string error = determineXMinMax(inputWS, xmins, xmaxs);
292 if (!error.empty())
293 throw std::runtime_error(error);
294
295 bool common_limits = true;
296 {
297 double xmin_common = xmins[0];
298 double xmax_common = xmaxs[0];
299 for (size_t i = 1; i < xmins.size(); ++i) {
300 if (xmins[i] != xmin_common) {
301 common_limits = false;
302 break;
303 }
304 if (xmaxs[i] != xmax_common) {
305 common_limits = false;
306 break;
307 }
308 }
309 }
310 if (common_limits) {
311 g_log.debug() << "Common limits between all spectra\n";
312 } else {
313 g_log.debug() << "Does not have common limits between all spectra\n";
314 }
315
316 // start doing actual work
317 EventWorkspace_const_sptr inputEventWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWS);
318 if (inputEventWS != nullptr) {
319 if (m_preserveEvents) {
320 if (inPlace) {
321 g_log.debug() << "Rebinning event workspace in place\n";
322 } else {
323 g_log.debug() << "Rebinning event workspace out of place\n";
324 outputWS = inputWS->clone();
325 }
326 auto outputEventWS = std::dynamic_pointer_cast<EventWorkspace>(outputWS);
327
328 if (common_limits) {
329 // get the delta from the first since they are all the same
330 BinEdges xValues(0);
331 const double delta = this->determineBinning(xValues.mutableRawData(), xmins[0], xmaxs[0]);
332 g_log.debug() << "delta = " << delta << "\n";
333 outputEventWS->setAllX(xValues);
334 } else {
335 // initialize progress reporting.
336 Progress prog(this, 0.0, 1.0, numSpectra);
337
338 // do the rebinning
339 PARALLEL_FOR_IF(Kernel::threadSafe(*inputEventWS, *outputWS))
340 for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
342 BinEdges xValues(0);
343 const double delta = this->determineBinning(xValues.mutableRawData(), xmins[wkspIndex], xmaxs[wkspIndex]);
344 g_log.debug() << "delta[wkspindex=" << wkspIndex << "] = " << delta << " xmin=" << xmins[wkspIndex]
345 << " xmax=" << xmaxs[wkspIndex] << "\n";
346 outputEventWS->setHistogram(wkspIndex, xValues);
347 prog.report(name()); // Report progress
349 }
351 }
352 } // end if (m_preserveEvents)
353 else // event workspace -> matrix workspace
354 {
355 //--------- Different output, OR you're inplace but not preserving Events
356 g_log.information() << "Creating a Workspace2D from the EventWorkspace " << inputEventWS->getName() << ".\n";
357 outputWS = create<DataObjects::Workspace2D>(*inputWS, numSpectra, HistogramData::BinEdges(m_numBins + 1));
358
359 // Initialize progress reporting.
360 Progress prog(this, 0.0, 1.0, numSpectra);
361
362 // Go through all the histograms and set the data
363 PARALLEL_FOR_IF(Kernel::threadSafe(*inputEventWS, *outputWS))
364 for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
366
367 // Set the X axis for each output histogram
368 MantidVec xValues;
369 const double delta = this->determineBinning(xValues, xmins[wkspIndex], xmaxs[wkspIndex]);
370 g_log.debug() << "delta[wkspindex=" << wkspIndex << "] = " << delta << "\n";
371 outputWS->setBinEdges(wkspIndex, xValues);
372
373 // Get a const event list reference. inputEventWS->dataY() doesn't work.
374 const EventList &el = inputEventWS->getSpectrum(wkspIndex);
375 MantidVec y_data, e_data;
376 // The EventList takes care of histogramming.
377 el.generateHistogram(xValues, y_data, e_data);
378
379 // Copy the data over.
380 outputWS->mutableY(wkspIndex) = y_data;
381 outputWS->mutableE(wkspIndex) = e_data;
382
383 // Report progress
384 prog.report(name());
386 }
388
389 // Copy all the axes
390 for (int i = 1; i < inputWS->axes(); i++) {
391 outputWS->replaceAxis(i, std::unique_ptr<Axis>(inputWS->getAxis(i)->clone(outputWS.get())));
392 outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit();
393 }
394
395 // Copy the units over too.
396 for (int i = 0; i < outputWS->axes(); ++i) {
397 outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit();
398 }
399 outputWS->setYUnit(inputEventWS->YUnit());
400 outputWS->setYUnitLabel(inputEventWS->YUnitLabel());
401 }
402 // Assign it to the output workspace property
403 setProperty("OutputWorkspace", outputWS);
404 return;
405 } else // (inputeventWS != NULL)
406 {
407 // workspace2d ----------------------------------------------------------
408 if (!m_isHistogram) {
409 g_log.information() << "Rebin: Converting Data to Histogram.\n";
410 Mantid::API::Algorithm_sptr ChildAlg = createChildAlgorithm("ConvertToHistogram");
411 ChildAlg->initialize();
412 ChildAlg->setProperty("InputWorkspace", inputWS);
413 ChildAlg->execute();
414 inputWS = ChildAlg->getProperty("OutputWorkspace");
415 }
416
417 // make output Workspace the same type is the input, but with new length of
418 // signal array
419 outputWS = API::WorkspaceFactory::Instance().create(inputWS, numSpectra, m_numBins + 1, m_numBins);
420
421 // Copy over the 'vertical' axis
422 if (inputWS->axes() > 1)
423 outputWS->replaceAxis(1, std::unique_ptr<Axis>(inputWS->getAxis(1)->clone(outputWS.get())));
424
425 Progress prog(this, 0.0, 1.0, numSpectra);
426 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
427 for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
429 // get const references to input Workspace arrays (no copying)
430 // TODO: replace with HistogramX/Y/E when VectorHelper::rebin is updated
431 const MantidVec &XValues = inputWS->readX(wkspIndex);
432 const MantidVec &YValues = inputWS->readY(wkspIndex);
433 const MantidVec &YErrors = inputWS->readE(wkspIndex);
434
435 // get references to output workspace data (no copying)
436 // TODO: replace with HistogramX/Y/E when VectorHelper::rebin is updated
437 MantidVec &YValues_new = outputWS->dataY(wkspIndex);
438 MantidVec &YErrors_new = outputWS->dataE(wkspIndex);
439
440 // create new output X axis
441 MantidVec XValues_new;
442 const double delta = this->determineBinning(XValues_new, xmins[wkspIndex], xmaxs[wkspIndex]);
443 g_log.debug() << "delta[wkspindex=" << wkspIndex << "] = " << delta << "\n";
444
445 // output data arrays are implicitly filled by function
446 try {
447 VectorHelper::rebin(XValues, YValues, YErrors, XValues_new, YValues_new, YErrors_new, m_isDistribution);
448 } catch (std::exception &ex) {
449 g_log.error() << "Error in rebin function: " << ex.what() << '\n';
450 throw;
451 }
452
453 // Populate the output workspace X values
454 outputWS->setBinEdges(wkspIndex, XValues_new);
455
456 prog.report(name());
458 }
460 outputWS->setDistribution(m_isDistribution);
461
462 // Now propagate any masking correctly to the output workspace
463 // More efficient to have this in a separate loop because
464 // MatrixWorkspace::maskBins blocks multi-threading
465 for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
466 if (inputWS->hasMaskedBins(wkspIndex)) // Does the current spectrum have any masked bins?
467 {
468 this->propagateMasks(inputWS, outputWS, wkspIndex);
469 }
470 }
471 // Copy the units over too.
472 for (int i = 0; i < outputWS->axes(); ++i) {
473 outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit();
474 }
475
476 if (!m_isHistogram) {
477 g_log.information() << "Rebin: Converting Data back to Data Points.\n";
478 Mantid::API::Algorithm_sptr ChildAlg = createChildAlgorithm("ConvertToPointData");
479 ChildAlg->initialize();
480 ChildAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", outputWS);
481 ChildAlg->execute();
482 outputWS = ChildAlg->getProperty("OutputWorkspace");
483 }
484
485 // Assign it to the output workspace property
486 setProperty("OutputWorkspace", outputWS);
487 } // end if (inputeventWS != NULL)
488}
489
490} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
double error
Definition: IndexPeaks.cpp:133
#define fabs(x)
Definition: Matrix.cpp:22
#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
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
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
ResampleX : TODO: DESCRIPTION.
Definition: ResampleX.h:20
void setOptions(const int numBins, const bool useLogBins, const bool isDist)
MADE PUBLIC FOR TESTING ONLY - DO NOT USE.
Definition: ResampleX.cpp:179
void exec() override
Execute the algorithm.
Definition: ResampleX.cpp:274
const std::string name() const override
Algorithm's name for identification.
Definition: ResampleX.cpp:34
const std::string category() const override
Algorithm's category for identification.
Definition: ResampleX.cpp:40
std::map< std::string, std::string > validateInputs() override
More complicated checks of parameters and their relations.
Definition: ResampleX.cpp:74
const std::string alias() const override
Algorithm's aliases.
Definition: ResampleX.cpp:42
double determineBinning(MantidVec &xValues, const double xmin, const double xmax)
MADE PUBLIC FOR TESTING ONLY - DO NOT USE.
Definition: ResampleX.cpp:194
void init() override
Initialize the algorithm's properties.
Definition: ResampleX.cpp:47
int version() const override
Algorithm's version for identification.
Definition: ResampleX.cpp:37
A class for holding :
Definition: EventList.h:56
void generateHistogram(const MantidVec &X, MantidVec &Y, MantidVec &E, bool skipError=false) const override
Generates both the Y and E (error) histograms w.r.t TOF for an EventList with or without WeightedEven...
Definition: EventList.cpp:2010
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.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
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
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Kernel::Logger g_log("ExperimentInfo")
static logger object
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
string determineXMinMax(const MatrixWorkspace_sptr &inputWS, vector< double > &xmins, vector< double > &xmaxs)
Determine the min and max x-values for each spectrum and error check the pairs.
Definition: ResampleX.cpp:110
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
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.
void MANTID_KERNEL_DLL rebin(const std::vector< double > &xold, const std::vector< double > &yold, const std::vector< double > &eold, const std::vector< double > &xnew, std::vector< double > &ynew, std::vector< double > &enew, bool distribution, bool addition=false)
Rebins data according to a new output X array.
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
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54