Newer
Older
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
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/InterpolatingRebin.h"
#include "MantidKernel/VectorHelper.h"
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
namespace Mantid
{
namespace Algorithms
{
// Register the class into the algorithm factory
DECLARE_ALGORITHM(InterpolatingRebin)
using namespace Kernel;
using namespace API;
/** Only calls its parent's (SimpleRebin) init()
*
*/
void InterpolatingRebin::init()
{
SimpleRebin::init();
}
/** Executes the rebin algorithm
*
* @throw runtime_error Thrown if the bin range does not intersect the range of the input workspace
*/
void InterpolatingRebin::exec()
{
// retrieve the properties
std::vector<double> rb_params=getProperty("Params");
// Get the input workspace
MatrixWorkspace_sptr inputW = getProperty("InputWorkspace");
//this calculation requires a distribution workspace but deal with the situation when we don't get this
const bool distCon = ! inputW->isDistribution();
if (distCon)
{
g_log.debug() << "Converting the input workspace to a distribution\n";
WorkspaceHelpers::makeDistribution(inputW);
}
DataObjects::Histogram1D::RCtype XValues_new;
// create new output X axis
const int ntcnew =
VectorHelper::createAxisFromRebinParams(rb_params,XValues_new.access());
const int nHists = inputW->getNumberHistograms();
// make output Workspace the same type is the input, but with new length of signal array
MatrixWorkspace_sptr outputW =
WorkspaceFactory::Instance().create(inputW, nHists, ntcnew, ntcnew-1);
// Copy over the 'vertical' axis
if (inputW->axes() > 1) outputW->replaceAxis( 1, inputW->getAxis(1)->clone(outputW.get()) );
//evaluate the rebinned data
outputXandEValues(inputW, XValues_new, outputW);
//check if there was a convert to distribution done previously
if (distCon)
{
g_log.debug() << "Converting the input and output workspaces _from_ distributions\n";
WorkspaceHelpers::makeDistribution(inputW, false);
// the calculation produces a distribution workspace but if they passed a non-distribution workspace they maybe not expect it, so convert back to the same form that was given
WorkspaceHelpers::makeDistribution(outputW, false);
}
outputW->isDistribution( ! distCon );
// Now propagate any masking correctly to the output workspace
// More efficient to have this in a separate loop because
// MatrixWorkspace::maskBins blocks multi-threading
for (int i=0; i < nHists; ++i)
{
if ( inputW->hasMaskedBins(i) ) // Does the current spectrum have any masked bins?
{
this->propagateMasks(inputW,outputW,i);
}
}
for (int i=0; i < outputW->axes(); ++i)
{
outputW->getAxis(i)->unit() = inputW->getAxis(i)->unit();
}
// Assign it to the output workspace property
setProperty("OutputWorkspace",outputW);
}
/** Calls the interpolation function for each histogram in the workspace
* @param[in] inputW workspace with un-interpolated data
* @param[in] XValues_new new x-values to interpolated to
* @param[out] outputW this will contain the interpolated data, the lengths of the histograms must corrospond with the number of x-values in XValues_new
*/
void InterpolatingRebin::outputXandEValues(API::MatrixWorkspace_const_sptr inputW, const DataObjects::Histogram1D::RCtype &XValues_new, API::MatrixWorkspace_sptr outputW)
{
g_log.debug() << "Preparing to calculate y-values using splines and estimate errors\n";
// prepare to use GSL functions but don't let them terminate Mantid
gsl_error_handler_t * old_handler = gsl_set_error_handler(NULL);
const int histnumber = inputW->getNumberHistograms();
Progress prog(this,0.0,1.0,histnumber);
PARALLEL_FOR2(inputW,outputW)
for (int hist=0; hist < histnumber;++hist)
{
PARALLEL_START_INTERUPT_REGION
// get const references to input Workspace arrays (no copying)
const MantidVec& XValues = inputW->readX(hist);
const MantidVec& YValues = inputW->readY(hist);
const MantidVec& YErrors = inputW->readE(hist);
//get references to output workspace data (no copying)
MantidVec& YValues_new=outputW->dataY(hist);
MantidVec& YErrors_new=outputW->dataE(hist);
try
{
// output data arrays are implicitly filled by function
cubicInterpolation(XValues, YValues, YErrors,
*XValues_new, YValues_new, YErrors_new);
}
catch (std::exception& ex)
{
g_log.error() << "Error in rebin function: " << ex.what() << std::endl;
throw;
}
// Populate the output workspace X values
outputW->setX(hist, XValues_new);
prog.report();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
gsl_set_error_handler(old_handler);
}
/**Uses cubic splines to interpolate the mean rate of change of the integral
* over the inputed data bins to that for the user supplied bins.
* Note that this algorithm was inplemented to provide a little more resolution
* on high count rate data. Whether it is more accurate than the standard rebin
* for all, or your, application needs more thought.
* The input data must be a distribution (proportional to the rate of change e.g.
* raw_counts/bin_widths) but note that these mean rate of counts data
* are integrals not (instanteously) sampled data. The error values on each point
* are a weighted mean of the error values from the surrounding input data. This
* makes sense if the interpolation error is low compared to the statistical
* errors on each input data point. The weighting is inversely proportional to
* the distance from the original data point to the new interpolated one.
*
* @param[in] xOld the x-values of the data that will be intepolated
* @param[in] yOld the data's y-values that corrospond to the x-values, must be 1 element shorter than xOld.
* @param[in] eOld the error on each y-value, must be same length as yOld.
* @param[in] xNew x-values to rebin to, must be monotonically increasing
* @param[out] yNew is overwritten with the algorithm output. Must be allocated and 1 element shorter than xnew.
* @param[out] eNew is overwritten with errors from the errors on the nearest input data points. Must be allocated with the same number of points as ynew
* @throw runtime_error if there is a problem executing one of the GSL functions
* @throw invalid_argument if any output x-values are outside the range of input x-values
**/
void InterpolatingRebin::cubicInterpolation(const MantidVec &xOld, const MantidVec &yOld, const MantidVec &eOld,
const MantidVec& xNew, MantidVec &yNew, MantidVec &eNew) const
{
// Make sure y and e vectors are of correct sizes
const size_t size_old = yOld.size();
if (size_old != (xOld.size() - 1) || size_old != eOld.size() )
throw std::runtime_error("rebin: y and error vectors should be of same size & 1 shorter than x");
const size_t size_new = yNew.size();
if (size_new != (xNew.size() - 1) || size_new != eNew.size() )
throw std::runtime_error("rebin: y and error vectors should be of same size & 1 shorter than x");
// get the bin centres of the input data
std::vector<double> xCensOld(size_new);
VectorHelper::convertToBinCentre(xOld, xCensOld);
// the centres of the output data
std::vector<double> xCensNew(size_new);
VectorHelper::convertToBinCentre(xNew, xCensNew);
// find the range of input values whose x-values just suround the output x-values
size_t oldIn1 =
std::lower_bound(xCensOld.begin(), xCensOld.end(), xCensNew.front())
- xCensOld.begin() - 1;
size_t oldIn2 =
std::lower_bound(xCensOld.begin(), xCensOld.end(), xCensNew.back())
- xCensOld.begin();
//bring one point before and one point after into the inpolation to reduce any errors coming in from the edge
oldIn1 --;
oldIn2 ++;
//check that the end points are all within the arrays
if ( oldIn1<0 || oldIn2>=size_old || oldIn1>oldIn2 )
{
throw std::invalid_argument("Problem with the requested x-values to intepolate to: There must be at\nleast two input data points below the range of intepolation points and\ntwo higher. Also the intepolation points must have monatomically increasing x-values.");
}
//get the GSL to allocate the memory, if this wasn't already done
gsl_interp_accel *acc = gsl_interp_accel_alloc();
const size_t nPoints = oldIn2 - oldIn1 + 1;
gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline, nPoints);
if ( ! acc || ! spline ||
//GSL calculates the splines
gsl_spline_init(spline, &xCensOld[oldIn1], &yOld[oldIn1], nPoints) )
{
throw std::runtime_error("Error setting up GSL spline functions");
}
for ( size_t i = 0; i < size_new; ++i )
{
yNew[i] = gsl_spline_eval(spline, xCensNew[i], acc);
//(basic) error estimate the based on a weighted mean of the errors of the surrounding input data points
eNew[i] = estimateError(xCensOld, eOld, xCensNew[i]);
}
//for GSL to clear up its memory use
gsl_spline_free (spline);
gsl_interp_accel_free (acc);
}
/**Estimates the error on each interpolated point by assuming it is similar to the errors in
* near by input data points. Output points with the same x-value as an input point have the
* same error as the input point. Points between input points have a error value that is a
* weighted mean of the closest input points
* @param[in] xsOld x-values of the input data around the point of interested
* @param[in] esOld error values for the same points in the input data as xsOld
* @param[in] xNew the value of x for at the point of interest
* @return the estimated error at that point
*/
double InterpolatingRebin::estimateError(const std::vector<double>& xsOld, const std::vector<double>& esOld,
const double xNew) const
{
//get the index of the first point that is higher in x, we'll base some of the error estimate on the error on this point
const int indAbove =
std::lower_bound(xsOld.begin(), xsOld.end(), xNew) - xsOld.begin();
const double error1 = esOld[indAbove];
// ratio of weightings will be inversely proportional to the distance between the points
double weight1 = xsOld[indAbove] - xNew;
//check if the points are close enough agnoring any spurious effects that can occur with exact comparisons of floating point numbers
if ( weight1 < 1e-100 )
{
// the point is on an input point, all the weight is on this point ignore the other
return error1;
}
weight1 = 1/weight1;
// if p were zero lower_bound must have found xCensNew <= xCensOld.front() but in that situation we should have exited before now
const double error2 = esOld[indAbove-1];
double weight2 = xNew - xsOld[indAbove-1];
if ( weight2 < 1e-100 )
{
// the point is on an input point, all the weight is on this point ignore the other
return error2;
}
weight2 = 1/weight2;
return ( weight1*error1 + weight2*error2 )/( weight1 + weight2 );
}
} // namespace Algorithm
} // namespace Mantid