/
SmoothData.cpp
197 lines (182 loc) · 7 KB
/
SmoothData.cpp
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
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/SmoothData.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/GroupingWorkspace.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/ArrayBoundedValidator.h"
namespace Mantid {
namespace Algorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(SmoothData)
using namespace Kernel;
using namespace API;
void SmoothData::init() {
declareProperty(
make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
"Name of the input workspace");
declareProperty(
make_unique<WorkspaceProperty<>>("OutputWorkspace", "",
Direction::Output),
"The name of the workspace to be created as the output of the algorithm");
std::vector<int> npts0;
npts0.push_back(3);
auto min = boost::make_shared<Kernel::ArrayBoundedValidator<int>>();
min->setLower(3);
// The number of points to use in the smoothing.
declareProperty(
Kernel::make_unique<ArrayProperty<int>>("NPoints", npts0, min,
Direction::Input),
"The number of points to average over (minimum 3). If an even number is\n"
"given, it will be incremented by 1 to make it odd (default value 3)");
declareProperty(
make_unique<WorkspaceProperty<Mantid::DataObjects::GroupingWorkspace>>(
"GroupingWorkspace", "", Direction::Input, PropertyMode::Optional),
"Optional: GroupingWorkspace to use for vector of NPoints.");
}
void SmoothData::exec() {
// Get the input properties
inputWorkspace = getProperty("InputWorkspace");
std::vector<int> nptsGroup = getProperty("NPoints");
Mantid::DataObjects::GroupingWorkspace_sptr groupWS =
getProperty("GroupingWorkspace");
if (groupWS) {
udet2group.clear();
int64_t nGroups;
groupWS->makeDetectorIDToGroupVector(udet2group, nGroups);
}
// Check that the number of points in the smoothing isn't larger than the
// spectrum length
const int vecSize = static_cast<int>(inputWorkspace->blocksize());
// Create the output workspace
MatrixWorkspace_sptr outputWorkspace =
WorkspaceFactory::Instance().create(inputWorkspace);
Progress progress(this, 0.0, 1.0, inputWorkspace->getNumberHistograms());
PARALLEL_FOR2(inputWorkspace, outputWorkspace)
// Loop over all the spectra in the workspace
for (int i = 0; i < static_cast<int>(inputWorkspace->getNumberHistograms());
++i) {
PARALLEL_START_INTERUPT_REGION
int npts = nptsGroup[0];
if (groupWS) {
const int group = validateSpectrumInGroup(static_cast<size_t>(i));
if (group < 0)
npts = 3;
else
npts = nptsGroup[group - 1];
}
if (npts >= vecSize) {
g_log.error("The number of averaging points requested is larger than the "
"spectrum length");
throw std::out_of_range("The number of averaging points requested is "
"larger than the spectrum length");
}
// Number of smoothing points must always be an odd number, so add 1 if it
// isn't.
if (!(npts % 2)) {
g_log.information("Adding 1 to number of smoothing points, since it must "
"always be odd");
++npts;
}
int halfWidth = (npts - 1) / 2;
// Copy the X data over. Preserves data sharing if present in input
// workspace.
outputWorkspace->setX(i, inputWorkspace->refX(i));
// Now get references to the Y & E vectors in the input and output
// workspaces
const MantidVec &Y = inputWorkspace->readY(i);
const MantidVec &E = inputWorkspace->readE(i);
MantidVec &newY = outputWorkspace->dataY(i);
MantidVec &newE = outputWorkspace->dataE(i);
if (npts == 0) {
newY = Y;
newE = E;
continue;
}
// Use total to help hold our moving average
double total = 0.0, totalE = 0.0;
// First push the values ahead of the current point onto total
for (int k = 0; k < halfWidth; ++k) {
if (Y[k] == Y[k])
total += Y[k]; // Exclude if NaN
totalE += E[k] * E[k];
}
// Now calculate the smoothed values for the 'end' points, where the number
// contributing
// to the smoothing will be less than NPoints
for (int j = 0; j <= halfWidth; ++j) {
const int index = j + halfWidth;
if (Y[index] == Y[index])
total += Y[index]; // Exclude if NaN
newY[j] = total / (index + 1);
totalE += E[index] * E[index];
newE[j] = sqrt(totalE) / (index + 1);
}
// This is the main part, where each data point is the average of NPoints
// points centred on the
// current point. Note that the statistical error will be reduced by
// sqrt(npts) because more
// data is now contributing to each point.
for (int k = halfWidth + 1; k < vecSize - halfWidth; ++k) {
const int kp = k + halfWidth;
const int km = k - halfWidth - 1;
total += (Y[kp] != Y[kp] ? 0.0 : Y[kp]) -
(Y[km] != Y[km] ? 0.0 : Y[km]); // Exclude if NaN
newY[k] = total / npts;
totalE += E[kp] * E[kp] - E[km] * E[km];
// Use of a moving average can lead to rounding error where what should be
// zero actually comes out as a tiny negative number - bad news for sqrt
// so protect
newE[k] = std::sqrt(std::abs(totalE)) / npts;
}
// This deals with the 'end' at the tail of each spectrum
for (int l = vecSize - halfWidth; l < vecSize; ++l) {
const int index = l - halfWidth;
total -=
(Y[index - 1] != Y[index - 1] ? 0.0 : Y[index - 1]); // Exclude if NaN
newY[l] = total / (vecSize - index);
totalE -= E[index - 1] * E[index - 1];
newE[l] = std::sqrt(std::abs(totalE)) / (vecSize - index);
}
progress.report();
PARALLEL_END_INTERUPT_REGION
} // Loop over spectra
PARALLEL_CHECK_INTERUPT_REGION
// Set the output workspace to its property
setProperty("OutputWorkspace", outputWorkspace);
}
//=============================================================================
/** Verify that all the contributing detectors to a spectrum belongs to the same
* group
* @param wi :: The workspace index in the workspace
* @return Group number if successful otherwise return -1
*/
int SmoothData::validateSpectrumInGroup(size_t wi) {
const std::set<detid_t> &dets =
inputWorkspace->getSpectrum(wi)->getDetectorIDs();
if (dets.empty()) // Not in group
{
g_log.debug() << wi << " <- this workspace index is empty!\n";
return -1;
}
auto it = dets.cbegin();
if (*it < 0) // bad pixel id
return -1;
try { // what if index out of range?
const int group = udet2group.at(*it);
if (group <= 0)
return -1;
++it;
for (; it != dets.end(); ++it) // Loop other all other udets
{
if (udet2group.at(*it) != group)
return -1;
}
return group;
} catch (...) {
}
return -1;
}
} // namespace Algorithms
} // namespace Mantid