-
Notifications
You must be signed in to change notification settings - Fork 122
/
MergeMD.cpp
328 lines (282 loc) · 11.8 KB
/
MergeMD.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
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
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidMDAlgorithms/MergeMD.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidDataObjects/MDBoxIterator.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/CPUTimer.h"
#include "MantidKernel/MandatoryValidator.h"
#include "MantidKernel/Strings.h"
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using namespace Mantid::DataObjects;
namespace Mantid {
namespace MDAlgorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MergeMD)
/// Algorithm's name for identification. @see Algorithm::name
const std::string MergeMD::name() const { return "MergeMD"; }
/// Algorithm's version for identification. @see Algorithm::version
int MergeMD::version() const { return 1; }
/// Algorithm's category for identification. @see Algorithm::category
const std::string MergeMD::category() const { return "MDAlgorithms\\Creation"; }
/** Initialize the algorithm's properties.
*/
void MergeMD::init() {
// declare arbitrary number of input m_workspaces as a list of strings at the
// moment
declareProperty(
std::make_unique<ArrayProperty<std::string>>(
"InputWorkspaces",
std::make_shared<MandatoryValidator<std::vector<std::string>>>()),
"The names of the input MDWorkspaces as a comma-separated list");
declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>(
"OutputWorkspace", "", Direction::Output),
"Name of the output MDWorkspace.");
// Set the box controller properties
this->initBoxControllerProps("2", 500, 16);
}
/** Create the output MDWorkspace from a list of input
*
* @param inputs :: list of names of input MDWorkspaces
*/
void MergeMD::createOutputWorkspace(std::vector<std::string> &inputs) {
auto it = inputs.begin();
for (; it != inputs.end(); it++) {
IMDEventWorkspace_sptr ws = std::dynamic_pointer_cast<IMDEventWorkspace>(
AnalysisDataService::Instance().retrieve(*it));
if (!ws)
throw std::invalid_argument(
"Workspace " + *it +
" is not a MDEventWorkspace. Cannot merge this workspace.");
else
m_workspaces.emplace_back(ws);
}
if (m_workspaces.empty())
throw std::invalid_argument("No valid m_workspaces specified.");
// Number of dimensions, start with the 0th one. They all have to match # of
// dims anyway
IMDEventWorkspace_const_sptr ws0 = m_workspaces[0];
size_t numDims = ws0->getNumDims();
auto displNorm = ws0->displayNormalization();
auto displNormH = ws0->displayNormalizationHisto();
// Extents to create.
std::vector<coord_t> dimMin(numDims, +1e30f);
std::vector<coord_t> dimMax(numDims, -1e30f);
// Validate each workspace
for (auto &ws : m_workspaces) {
if (ws->getNumDims() != numDims)
throw std::invalid_argument(
"Workspace " + ws->getName() +
" does not match the number of dimensions of the others (" +
Strings::toString(ws->getNumDims()) + ", expected " +
Strings::toString(numDims) + ")");
if (ws->getEventTypeName() != ws0->getEventTypeName())
throw std::invalid_argument(
"Workspace " + ws->getName() +
" does not match the MDEvent type of the others (" +
ws->getEventTypeName() + ", expected " + ws0->getEventTypeName() +
")");
for (size_t d = 0; d < numDims; d++) {
IMDDimension_const_sptr dim = ws->getDimension(d);
IMDDimension_const_sptr dim0 = ws0->getDimension(d);
if (dim->getName() != dim0->getName())
throw std::invalid_argument("Workspace " + ws->getName() +
" does not have the same dimension " +
Strings::toString(d) + " as the others (" +
dim->getName() + ", expected " +
dim0->getName() + ")");
// Find the extents
if (dim->getMaximum() > dimMax[d])
dimMax[d] = dim->getMaximum();
if (dim->getMinimum() < dimMin[d])
dimMin[d] = dim->getMinimum();
}
}
// OK, now create the blank MDWorkspace
// Have the factory create it
out = MDEventFactory::CreateMDWorkspace(numDims, ws0->getEventTypeName());
out->setDisplayNormalization(displNorm);
out->setDisplayNormalizationHisto(displNormH);
// Give all the dimensions
for (size_t d = 0; d < numDims; d++) {
IMDDimension_const_sptr dim0 = ws0->getDimension(d);
MDHistoDimension *dim = new MDHistoDimension(
dim0->getName(), dim0->getDimensionId(), dim0->getMDFrame(), dimMin[d],
dimMax[d], dim0->getNBins());
out->addDimension(MDHistoDimension_sptr(dim));
}
// Initialize it using the dimension
out->initialize();
// Set the box controller settings from the properties
this->setBoxController(out->getBoxController());
// Perform the initial box splitting
out->splitBox();
// copy experiment infos
uint16_t nExperiments(0);
if (m_workspaces.size() > std::numeric_limits<uint16_t>::max())
throw std::invalid_argument(
"currently we can not combine more then 65535 experiments");
else
nExperiments = static_cast<uint16_t>(m_workspaces.size());
for (uint16_t i = 0; i < nExperiments; i++) {
uint16_t nWSexperiments = m_workspaces[i]->getNumExperimentInfo();
experimentInfoNo.emplace_back(nWSexperiments);
for (uint16_t j = 0; j < nWSexperiments; j++) {
API::ExperimentInfo_sptr ei = API::ExperimentInfo_sptr(
m_workspaces[i]->getExperimentInfo(j)->cloneExperimentInfo());
out->addExperimentInfo(ei);
}
}
// Cumulative sum of number of experimentInfo and reverse order
std::partial_sum(experimentInfoNo.begin(), experimentInfoNo.end(),
experimentInfoNo.begin());
std::reverse(std::begin(experimentInfoNo), std::end(experimentInfoNo));
}
//----------------------------------------------------------------------------------------------
/** Copy the extra data (not signal, error or coordinates) from one event to
* another with different numbers of dimensions
*
* @param srcEvent :: the source event, being copied
* @param newEvent :: the destination event
* @param runIndexOffset :: offset to be added to the runIndex
*/
template <size_t nd, size_t ond>
inline void copyEvent(const MDLeanEvent<nd> &srcEvent,
MDLeanEvent<ond> &newEvent,
const uint16_t runIndexOffset) {
// Nothing extra copy - this is no-op
UNUSED_ARG(srcEvent);
UNUSED_ARG(newEvent);
UNUSED_ARG(runIndexOffset);
}
//----------------------------------------------------------------------------------------------
/** Copy the extra data (not signal, error or coordinates) from one event to
* another with different numbers of dimensions
*
* @param srcEvent :: the source event, being copied
* @param newEvent :: the destination event
* @param runIndexOffset :: offset to be added to the runIndex
*/
template <size_t nd, size_t ond>
inline void copyEvent(const MDEvent<nd> &srcEvent, MDEvent<ond> &newEvent,
const uint16_t runIndexOffset) {
newEvent.setDetectorId(srcEvent.getDetectorID());
newEvent.setRunIndex(
static_cast<uint16_t>(srcEvent.getRunIndex() + runIndexOffset));
}
//----------------------------------------------------------------------------------------------
/** Perform the adding.
* Will do out += ws
*
* @param ws :: MDEventWorkspace to clone
*/
template <typename MDE, size_t nd>
void MergeMD::doPlus(typename MDEventWorkspace<MDE, nd>::sptr ws2) {
// CPUTimer tim;
typename MDEventWorkspace<MDE, nd>::sptr ws1 =
std::dynamic_pointer_cast<MDEventWorkspace<MDE, nd>>(out);
if (!ws1 || !ws2)
throw std::runtime_error("Incompatible workspace types passed to MergeMD.");
MDBoxBase<MDE, nd> *box1 = ws1->getBox();
MDBoxBase<MDE, nd> *box2 = ws2->getBox();
uint16_t runIndexOffset = experimentInfoNo.back();
experimentInfoNo.pop_back();
// How many events you started with
size_t initial_numEvents = ws1->getNPoints();
// Make a leaf-only iterator through all boxes with events in the RHS
// workspace
std::vector<API::IMDNode *> boxes;
box2->getBoxes(boxes, 1000, true);
auto numBoxes = int(boxes.size());
bool fileBasedSource(false);
if (ws2->isFileBacked())
fileBasedSource = true;
// Add the boxes in parallel. They should be spread out enough on each
// core to avoid stepping on each other.
// cppcheck-suppress syntaxError
PRAGMA_OMP( parallel for if (!ws2->isFileBacked()) )
for (int i = 0; i < numBoxes; i++) {
PARALLEL_START_INTERUPT_REGION
auto *box = dynamic_cast<MDBox<MDE, nd> *>(boxes[i]);
if (box && !box->getIsMasked()) {
// Copy the events from WS2 and add them into WS1
const std::vector<MDE> &events = box->getConstEvents();
// Add events, with bounds checking
for (auto it = events.cbegin(); it != events.cend(); ++it) {
// Create the event
MDE newEvent(it->getSignal(), it->getErrorSquared(), it->getCenter());
// Copy extra data, if any
copyEvent(*it, newEvent, runIndexOffset);
// Add it to the workspace
box1->addEvent(newEvent);
}
if (fileBasedSource)
box->clear();
else
box->releaseEvents();
}
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
// Progress * prog2 = new Progress(this, 0.4, 0.9, 100);
Progress *prog2 = nullptr;
ThreadScheduler *ts = new ThreadSchedulerFIFO();
ThreadPool tp(ts, 0, prog2);
ws1->splitAllIfNeeded(ts);
// prog2->resetNumSteps( ts->size(), 0.4, 0.6);
tp.joinAll();
// Set a marker that the file-back-end needs updating if the # of events
// changed.
if (ws1->getNPoints() != initial_numEvents)
ws1->setFileNeedsUpdating(true);
//
// std::cout << tim << " to add workspace " << ws2->name() << '\n';
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void MergeMD::exec() {
CPUTimer tim;
// Check that all input workspaces exist and match in certain important ways
const std::vector<std::string> inputs_orig = getProperty("InputWorkspaces");
// This will hold the inputs, with the groups separated off
std::vector<std::string> inputs;
for (const auto &input : inputs_orig) {
WorkspaceGroup_sptr wsgroup =
AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>(input);
if (wsgroup) { // Workspace group
std::vector<std::string> group = wsgroup->getNames();
inputs.insert(inputs.end(), group.begin(), group.end());
} else {
// Single workspace
inputs.emplace_back(input);
}
}
if (inputs.size() == 1) {
g_log.error("Only one input workspace specified");
throw std::invalid_argument("Only one input workspace specified");
}
// Create a blank output workspace
this->createOutputWorkspace(inputs);
// Run PlusMD on each of the input workspaces, in order.
double progStep = 1.0 / double(m_workspaces.size());
for (size_t i = 0; i < m_workspaces.size(); i++) {
g_log.information() << "Adding workspace " << m_workspaces[i]->getName()
<< '\n';
progress(double(i) * progStep, m_workspaces[i]->getName());
CALL_MDEVENT_FUNCTION(doPlus, m_workspaces[i]);
}
this->progress(0.95, "Refreshing cache");
out->refreshCache();
this->setProperty("OutputWorkspace", out);
g_log.debug() << tim << " to merge all workspaces.\n";
}
} // namespace MDAlgorithms
} // namespace Mantid