Mantid
Loading...
Searching...
No Matches
InterpolatingRebin.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"
13#include "MantidHistogramData/Histogram.h"
17
18#include <gsl/gsl_errno.h>
19#include <gsl/gsl_interp.h>
20#include <gsl/gsl_spline.h>
21
22#include <boost/lexical_cast.hpp>
23
24namespace Mantid::Algorithms {
25
26// Register the class into the algorithm factory
27DECLARE_ALGORITHM(InterpolatingRebin)
28
29using namespace Kernel;
30using namespace API;
31using namespace HistogramData;
32using namespace DataObjects;
33
38 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
39 "Workspace containing the input data");
40 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
41 "The name to give the output workspace");
42
43 declareProperty(std::make_unique<ArrayProperty<double>>("Params", std::make_shared<RebinParamsValidator>()),
44 "A comma separated list of first bin boundary, width, last bin boundary. "
45 "Optionally "
46 "this can be followed by a comma and more widths and last boundary "
47 "pairs. "
48 "Optionally this can also be a single number, which is the bin width. "
49 "In this case, the boundary of binning will be determined by minimum and "
50 "maximum TOF "
51 "values among all events, or previous binning boundary, in case of event "
52 "Workspace, or "
53 "non-event Workspace, respectively. Negative width values indicate "
54 "logarithmic binning. ");
55}
56
63 // Get the input workspace
64 MatrixWorkspace_sptr inputW = getProperty("InputWorkspace");
65
66 // retrieve the properties
67 std::vector<double> rb_params = Rebin::rebinParamsFromInput(getProperty("Params"), *inputW, g_log);
68 HistogramData::BinEdges XValues_new(0);
69 // create new output X axis
70 const int ntcnew = VectorHelper::createAxisFromRebinParams(rb_params, XValues_new.mutableRawData());
71
72 const auto nHists = static_cast<int>(inputW->getNumberHistograms());
73 // make output Workspace the same type as the input but with the new axes
74 MatrixWorkspace_sptr outputW = create<MatrixWorkspace>(*inputW, BinEdges(ntcnew));
75 // Copy over the 'vertical' axis
76 if (inputW->axes() > 1)
77 outputW->replaceAxis(1, std::unique_ptr<Axis>(inputW->getAxis(1)->clone(outputW.get())));
78 outputW->setDistribution(true);
79
80 // this calculation requires a distribution workspace but deal with the
81 // situation when we don't get this
82 bool distCon(false);
83 if (!inputW->isDistribution()) {
84 g_log.debug() << "Converting the input workspace to a distribution\n";
86 distCon = true;
87 }
88
89 try {
90 // evaluate the rebinned data
91 outputYandEValues(inputW, XValues_new, outputW);
92 } catch (std::exception &) {
93
94 if (distCon) {
95 // we need to return the input workspace to the state we found it in
97 }
98 throw;
99 }
100
101 // check if there was a convert to distribution done previously
102 if (distCon) {
103 g_log.debug() << "Converting the input and output workspaces _from_ distributions\n";
105 // the calculation produces a distribution workspace but if they passed a
106 // non-distribution workspace they maybe not expect it, so convert back to
107 // the same form that was given
109 outputW->setDistribution(false);
110 }
111
112 // Now propagate any masking correctly to the output workspace
113 // More efficient to have this in a separate loop because
114 // MatrixWorkspace::maskBins blocks multi-threading
115 for (int i = 0; i < nHists; ++i) {
116 if (inputW->hasMaskedBins(i)) // Does the current spectrum have any masked bins?
117 {
118 this->propagateMasks(inputW, outputW, i);
119 }
120 }
121
122 for (int i = 0; i < outputW->axes(); ++i) {
123 outputW->getAxis(i)->unit() = inputW->getAxis(i)->unit();
124 }
125
126 // Assign it to the output workspace property
127 setProperty("OutputWorkspace", outputW);
128}
136 const HistogramData::BinEdges &XValues_new,
137 const API::MatrixWorkspace_sptr &outputW) {
138 g_log.debug() << "Preparing to calculate y-values using splines and estimate errors\n";
139
140 // prepare to use GSL functions but don't let them terminate Mantid
141 gsl_error_handler_t *old_handler = gsl_set_error_handler(nullptr);
142
143 const auto histnumber = static_cast<int>(inputW->getNumberHistograms());
144 Progress prog(this, 0.0, 1.0, histnumber);
145 for (int hist = 0; hist < histnumber; ++hist) {
146
147 try {
148 // output data arrays are implicitly filled by function
149 outputW->setHistogram(hist, cubicInterpolation(inputW->histogram(hist), XValues_new));
150 } catch (std::exception &ex) {
151 g_log.error() << "Error in rebin function: " << ex.what() << '\n';
152 throw;
153 }
154
155 // Populate the output workspace X values
156 outputW->setBinEdges(hist, XValues_new);
157
158 prog.report();
159 }
160
161 gsl_set_error_handler(old_handler);
162}
163
190Histogram InterpolatingRebin::cubicInterpolation(const Histogram &oldHistogram, const BinEdges &xNew) const {
191 const auto &yOld = oldHistogram.y();
192
193 const size_t size_old = yOld.size();
194 if (size_old == 0)
195 throw std::runtime_error("Empty spectrum found, aborting!");
196
197 const size_t size_new = xNew.size() - 1; // -1 because BinEdges
198
199 // get the bin centres of the input data
200 auto xCensOld = oldHistogram.points();
201 VectorHelper::convertToBinCentre(oldHistogram.x().rawData(), xCensOld.mutableRawData());
202 // the centres of the output data
203 Points xCensNew(size_new);
204 VectorHelper::convertToBinCentre(xNew.rawData(), xCensNew.mutableRawData());
205
206 // find the range of input values whose x-values just suround the output
207 // x-values
208 size_t oldIn1 = std::lower_bound(xCensOld.begin(), xCensOld.end(), xCensNew.front()) - xCensOld.begin();
209 if (oldIn1 == 0) { // the lowest interpolation value might be out of range but
210 // if it is almost on the boundary let it through
211 if (std::abs(xCensOld.front() - xCensNew.front()) < 1e-8 * (xCensOld.back() - xCensOld.front())) {
212 oldIn1 = 1;
213 // make what should be a very small correction
214 xCensNew.mutableRawData().front() = xCensOld.front();
215 }
216 }
217
218 size_t oldIn2 = std::lower_bound(xCensOld.begin(), xCensOld.end(), xCensNew.back()) - xCensOld.begin();
219 if (oldIn2 == size_old) { // the highest point is nearly out of range of the
220 // input data but if it's very near the boundary let
221 // it through
222 if (std::abs(xCensOld.back() - xCensNew.back()) < 1e-8 * (xCensOld.back() - xCensOld.front())) {
223 oldIn2 = size_old - 1;
224 // make what should be a very small correction
225 xCensNew.mutableRawData().back() = xCensOld.back();
226 }
227 }
228
229 // check that the intepolation points fit well enough within the data for
230 // reliable intepolation to be done
231 bool goodRangeLow(false), goodRangeHigh(false), canInterpol(false);
232 if (oldIn1 > 1) { // set the range of the fit, including more input data to
233 // improve accuracy
234 oldIn1 -= 2;
235 goodRangeLow = true;
236 canInterpol = true;
237 } else {
238 if (oldIn1 > 0) {
239 canInterpol = true;
240 oldIn1--;
241 }
242 }
243
244 if (oldIn2 < size_old - 1) {
245 oldIn2++;
246 goodRangeHigh = true;
247 } else {
248 if (oldIn2 >= size_old) {
249 canInterpol = false;
250 }
251 }
252
253 const auto &xOld = oldHistogram.x();
254 const auto &eOld = oldHistogram.e();
255
256 // No Interpolation branch
257 if (!canInterpol) {
258 if (VectorHelper::isConstantValue(yOld.rawData())) {
259 // this copies the single y-value into the output array, errors are still
260 // calculated from the nearest input data points
261 // this is as much as we need to do in this (trival) case
262 return noInterpolation(oldHistogram, xNew);
263 } else { // some points are two close to the edge of the data
264
265 throw std::invalid_argument(std::string("At least one x-value to interpolate to is outside the "
266 "range of the original data.\n") +
267 "original data range: " + boost::lexical_cast<std::string>(xOld.front()) + " to " +
268 boost::lexical_cast<std::string>(xOld.back()) + "\n" +
269 "range to try to interpolate to " + boost::lexical_cast<std::string>(xNew.front()) +
270 " to " + boost::lexical_cast<std::string>(xNew.back()));
271 }
272 }
273
274 // Can Interpolate
275 // Create Histogram
276 Histogram newHistogram{xNew, Frequencies(xNew.size() - 1)};
277
278 auto &yNew = newHistogram.mutableY();
279 auto &eNew = newHistogram.mutableE();
280
281 if ((!goodRangeLow) || (!goodRangeHigh)) {
282 g_log.information() << "One or more points in the interpolation are near "
283 "the boundary of the input data, these points will "
284 "have slightly less accuracy\n";
285 }
286
287 // get the GSL to allocate the memory
288 gsl_interp_accel *acc = nullptr;
289 gsl_spline *spline = nullptr;
290 try {
291 acc = gsl_interp_accel_alloc();
292 const size_t nPoints = oldIn2 - oldIn1 + 1;
293 spline = gsl_spline_alloc(gsl_interp_cspline, nPoints);
294
295 // test the allocation
296 if (!acc || !spline ||
297 // calculate those splines, GSL uses pointers to the vector array (which
298 // is always contiguous)
299 gsl_spline_init(spline, &xCensOld[oldIn1], &yOld[oldIn1], nPoints)) {
300 throw std::runtime_error("Error setting up GSL spline functions");
301 }
302
303 for (size_t i = 0; i < size_new; ++i) {
304 yNew[i] = gsl_spline_eval(spline, xCensNew[i], acc);
305 //(basic) error estimate the based on a weighted mean of the errors of the
306 // surrounding input data points
307 eNew[i] = estimateError(xCensOld, eOld, xCensNew[i]);
308 }
309 }
310 // for GSL to clear up its memory use
311 catch (std::exception &) {
312 if (acc) {
313 if (spline) {
314 gsl_spline_free(spline);
315 }
316 gsl_interp_accel_free(acc);
317 }
318 throw;
319 }
320 gsl_spline_free(spline);
321 gsl_interp_accel_free(acc);
322
323 return newHistogram;
324}
334Histogram InterpolatingRebin::noInterpolation(const Histogram &oldHistogram, const BinEdges &xNew) const {
335 Histogram newHistogram{xNew, Frequencies(xNew.size() - 1)};
336 auto &yNew = newHistogram.mutableY();
337 auto &eNew = newHistogram.mutableE();
338
339 yNew.assign(yNew.size(), oldHistogram.y().front());
340
341 const auto &xPointData = oldHistogram.points();
342 const auto &eOld = oldHistogram.e();
343
344 // -1 because xNew.size is 1 bigger than eNew
345 std::transform(xNew.cbegin(), xNew.cend() - 1, eNew.begin(),
346 [&](double x) { return estimateError(xPointData, eOld, x); });
347
348 return newHistogram;
349}
361double InterpolatingRebin::estimateError(const Points &xsOld, const HistogramE &esOld, const double xNew) const {
362
363 // find the first point in the array that has a higher value of x, we'll base
364 // some of the error estimate on the error on this point
365
366 const size_t indAbove = std::lower_bound(xsOld.begin(), xsOld.end(), xNew) - xsOld.begin();
367
368 // if the point's x-value is out of the range covered by the x-values in the
369 // input data return the error value at the end of the range
370 if (indAbove == 0) {
371 return esOld.front();
372 }
373 // xsOld is 1 longer than xsOld
374 if (indAbove >= esOld.size()) { // cubicInterpolation() checks that that there
375 // are no empty histograms
376 return esOld.back();
377 }
378
379 const double error1 = esOld[indAbove];
380 // ratio of weightings will be inversely proportional to the distance between
381 // the points
382 double weight1 = xsOld[indAbove] - xNew;
383 // check if the points are close enough agnoring any spurious effects that can
384 // occur with exact comparisons of floating point numbers
385 if (weight1 < 1e-100) {
386 // the point is on an input point, all the weight is on this point ignore
387 // the other
388 return error1;
389 }
390 weight1 = 1 / weight1;
391
392 // if p were zero lower_bound must have found xCensNew <= xCensOld.front() but
393 // in that situation we should have exited before now
394 const double error2 = esOld[indAbove - 1];
395 double weight2 = xNew - xsOld[indAbove - 1];
396 if (weight2 < 1e-100) {
397 // the point is on an input point, all the weight is on this point ignore
398 // the other
399 return error2;
400 }
401 weight2 = 1 / weight2;
402
403 return (weight1 * error1 + weight2 * error2) / (weight1 + weight2);
404}
405
406} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void outputYandEValues(const API::MatrixWorkspace_const_sptr &inputW, const HistogramData::BinEdges &XValues_new, const API::MatrixWorkspace_sptr &outputW)
Calls the interpolation function for each histogram in the workspace.
void exec() override
Executes the rebin algorithm.
HistogramData::Histogram noInterpolation(const HistogramData::Histogram &oldHistogram, const HistogramData::BinEdges &xNew) const
This can be used whenever the original spectrum is filled with only one value.
double estimateError(const HistogramData::Points &xsOld, const HistogramData::HistogramE &esOld, const double xNew) const
Estimates the error on each interpolated point by assuming it is similar to the errors in near by inp...
HistogramData::Histogram cubicInterpolation(const HistogramData::Histogram &oldHistogram, const HistogramData::BinEdges &xNew) const
Uses cubic splines to interpolate the mean rate of change of the integral over the inputed data bins ...
void init() override
Only calls its parent's (Rebin) init()
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
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
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 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
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
void MANTID_KERNEL_DLL convertToBinCentre(const std::vector< double > &bin_edges, std::vector< double > &bin_centres)
Convert an array of bin boundaries to bin center values.
bool MANTID_KERNEL_DLL isConstantValue(const std::vector< double > &arra)
Assess if all the values in the vector are equal or if there are some different values.
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.
static void makeDistribution(const MatrixWorkspace_sptr &workspace, const bool forwards=true)
Divides the data in a workspace by the bin width to make it a distribution.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54