13#include "MantidHistogramData/Histogram.h"
18#include <gsl/gsl_errno.h>
19#include <gsl/gsl_interp.h>
20#include <gsl/gsl_spline.h>
22#include <boost/lexical_cast.hpp>
29using namespace Kernel;
31using namespace HistogramData;
32using namespace DataObjects;
39 "Workspace containing the input data");
41 "The name to give the output workspace");
44 "A comma separated list of first bin boundary, width, last bin boundary. "
46 "this can be followed by a comma and more widths and last boundary "
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 "
51 "values among all events, or previous binning boundary, in case of event "
53 "non-event Workspace, respectively. Negative width values indicate "
54 "logarithmic binning. ");
68 HistogramData::BinEdges XValues_new(0);
72 const auto nHists =
static_cast<int>(inputW->getNumberHistograms());
76 if (inputW->axes() > 1)
77 outputW->replaceAxis(1, std::unique_ptr<Axis>(inputW->getAxis(1)->clone(outputW.get())));
78 outputW->setDistribution(
true);
83 if (!inputW->isDistribution()) {
84 g_log.
debug() <<
"Converting the input workspace to a distribution\n";
92 }
catch (std::exception &) {
103 g_log.
debug() <<
"Converting the input and output workspaces _from_ distributions\n";
109 outputW->setDistribution(
false);
115 for (
int i = 0; i < nHists; ++i) {
116 if (inputW->hasMaskedBins(i))
122 for (
int i = 0; i < outputW->axes(); ++i) {
123 outputW->getAxis(i)->unit() = inputW->getAxis(i)->unit();
136 const HistogramData::BinEdges &XValues_new,
138 g_log.
debug() <<
"Preparing to calculate y-values using splines and estimate errors\n";
141 gsl_error_handler_t *old_handler = gsl_set_error_handler(
nullptr);
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) {
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';
156 outputW->setBinEdges(hist, XValues_new);
161 gsl_set_error_handler(old_handler);
191 const auto &yOld = oldHistogram.y();
193 const size_t size_old = yOld.size();
195 throw std::runtime_error(
"Empty spectrum found, aborting!");
197 const size_t size_new = xNew.size() - 1;
200 auto xCensOld = oldHistogram.points();
203 Points xCensNew(size_new);
208 size_t oldIn1 = std::lower_bound(xCensOld.begin(), xCensOld.end(), xCensNew.front()) - xCensOld.begin();
211 if (std::abs(xCensOld.front() - xCensNew.front()) < 1e-8 * (xCensOld.back() - xCensOld.front())) {
214 xCensNew.mutableRawData().front() = xCensOld.front();
218 size_t oldIn2 = std::lower_bound(xCensOld.begin(), xCensOld.end(), xCensNew.back()) - xCensOld.begin();
219 if (oldIn2 == size_old) {
222 if (std::abs(xCensOld.back() - xCensNew.back()) < 1e-8 * (xCensOld.back() - xCensOld.front())) {
223 oldIn2 = size_old - 1;
225 xCensNew.mutableRawData().back() = xCensOld.back();
231 bool goodRangeLow(
false), goodRangeHigh(
false), canInterpol(
false);
244 if (oldIn2 < size_old - 1) {
246 goodRangeHigh =
true;
248 if (oldIn2 >= size_old) {
253 const auto &xOld = oldHistogram.x();
254 const auto &eOld = oldHistogram.e();
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()));
276 Histogram newHistogram{xNew, Frequencies(xNew.size() - 1)};
278 auto &yNew = newHistogram.mutableY();
279 auto &eNew = newHistogram.mutableE();
281 if ((!goodRangeLow) || (!goodRangeHigh)) {
283 "the boundary of the input data, these points will "
284 "have slightly less accuracy\n";
288 gsl_interp_accel *acc =
nullptr;
289 gsl_spline *spline =
nullptr;
291 acc = gsl_interp_accel_alloc();
292 const size_t nPoints = oldIn2 - oldIn1 + 1;
293 spline = gsl_spline_alloc(gsl_interp_cspline, nPoints);
296 if (!acc || !spline ||
299 gsl_spline_init(spline, &xCensOld[oldIn1], &yOld[oldIn1], nPoints)) {
300 throw std::runtime_error(
"Error setting up GSL spline functions");
303 for (
size_t i = 0; i < size_new; ++i) {
304 yNew[i] = gsl_spline_eval(spline, xCensNew[i], acc);
311 catch (std::exception &) {
314 gsl_spline_free(spline);
316 gsl_interp_accel_free(acc);
320 gsl_spline_free(spline);
321 gsl_interp_accel_free(acc);
335 Histogram newHistogram{xNew, Frequencies(xNew.size() - 1)};
336 auto &yNew = newHistogram.mutableY();
337 auto &eNew = newHistogram.mutableE();
339 yNew.assign(yNew.size(), oldHistogram.y().front());
341 const auto &xPointData = oldHistogram.points();
342 const auto &eOld = oldHistogram.e();
345 std::transform(xNew.cbegin(), xNew.cend() - 1, eNew.begin(),
346 [&](
double x) { return estimateError(xPointData, eOld, x); });
366 const size_t indAbove = std::lower_bound(xsOld.begin(), xsOld.end(), xNew) - xsOld.begin();
371 return esOld.front();
374 if (indAbove >= esOld.size()) {
379 const double error1 = esOld[indAbove];
382 double weight1 = xsOld[indAbove] - xNew;
385 if (weight1 < 1e-100) {
390 weight1 = 1 / weight1;
394 const double error2 = esOld[indAbove - 1];
395 double weight2 = xNew - xsOld[indAbove - 1];
396 if (weight2 < 1e-100) {
401 weight2 = 1 / weight2;
403 return (weight1 * error1 + weight2 * error2) / (weight1 + weight2);
#define DECLARE_ALGORITHM(classname)
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Helper class for reporting progress from algorithms.
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...
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.
Support for a property that holds an array of values.
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.
void error(const std::string &msg)
Logs at error level.
void information(const std::string &msg)
Logs at information level.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
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 > ¶ms, 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.
@ Output
An output workspace.